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In-flight  simulations  are  normally  accomplished  by  using 
model-following  control  laws  which  depend  on  accurate  knowledge  of  the 
stability  derivatives  of  the  host  aircraft.  Degraded  simulation 
performance  may  result  if  the  stability  derivatives  deviate  considerably 
from  their  presumed  values.  Gain  scheduling  is  often  employed  to 
compensate  for  plant  parameter  variations,  but  this  form  of  open-loop 
compensation  usually  requires  extensive  flight  testing  for  proper  fine 
tuning.  This  thesis  implements  an  adaptive,  fast-sampling  control  law 
to  compensate  for  changing  aircraft  parameters.  The  step-response 
matrix  which  is  required  for  this  implementation  is  identified 
recursively  using  a  recently  developed  technique  which  does  not  require 
special  "test"  signals  and  which  automatically  discounts  old  data 
depending  on  the  input  excitation  detected.  Tracking  fidelity  is 
maintained  despite  parameter  changes  which  occur  either  abruptly  or 
slowly.  Simulations  are  conducted,  using  a  model  of  the  AFTI/F-16 
aircraft  and  the  control  design  package  MATRIX^,  to  test  the  resulting 
adaptive  system.  Actuator  position  and  rate  limits  are  discussed.  The 
performance  of  the  resulting  system  is  excellent  and  demonstrates  the 
relative  advantages  of  adaptive  controllers  for  in-flight  simulation. 
Recommendations  are  made  for  future  analysis  including  the  use  of 
moving-bank  estimators. 


xxxii 


PARAMETER-ADAPTIVE  MODEL-FOLLOWING  FOR  IN-FLIGHT  SIMULATION 


I.  Introduction 
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Flight  simulation  plays  an  essential  role  in  the  development  of 
modern  day  aircraft  as  well  as  pilot  training.  It  can  be  conducted  in  a 
ground  based  simulator,  or  in  actual  flight,  with  a  special  class  of 
aircraft  called  in-flight  simulators.  An  in-flight  simulator  is  an 
aircraft  whose  stability,  feel,  and  flying  characteristics  can  be 
changed  to  match  those  of  another  aircraft.  This  is  accomplished  by 
interfacing  the  pilot  to  the  aircraft  by  means  of  a  "fly-by-wire" 
variable  stability  flight  control  system  and  a  programmable  artificial 
feel  system.  As  the  pilot  moves  the  controls,  the  aircraft  responds  as 
the  simulated  vehicle  would.  The  pilot  experiences  the  real  flight 
motions  and  handling  qualities  of  the  simulated  aircraft  (29:38). 
In-flight  simulators  are  essential  tools  in  the  research  and  development 
process  that  provide  the  capability  to  realistically  and  safely  evaluate 
new  or  modified  aircraft  and  weapon  systems  before  first-flight  and 
before  committing  to  production  (9).  They  have  been  used  in  the 
development  of  the  F-16,  YF-17,  F-18,  A-10,  B-1,  and  Space  Shuttle. 
Other  uses  invole  many  types  of  generic  research,  especially  handling 
qualities,  and  for  several  types  of  specialized  training,  particularly 
at  the  Air  Force  and  Navy  Test  Pilot  Schools  (29:38). 


Simulator  Test  Aircraft),  will  meet  future  research  needs  by  providing 
upgraded  simulation  capabilities  (Figure  1-2).  This  will  be 
accomplished  by  modifying  a  modern,  high  performance  fighter  aircraft 
with  variable  stability  controls  and  a  reprogrammable  cockpit  (Figure 
1-3). 


1.2  Problem  Description 

The  dominant  performance  features  of  the  VISTA/F-16  are  likely 
to  be  a  wide  range  of  operating  conditions  in  an  expanded  flight 
envelope  (Figure  1-4),  and  large  amplitude  maneuvering  (8,42).  Both  of 
these  situations  produce  large  variations  in  the  stability  derivatives 
of  the  host  aircraft.  This  variation  of  parameters,  coupled  with 
initial  uncertainty  on  their  exact  values,  may  adversely  affect  the 
fidelity  of  in-flight  simulations  (given  a  fixed  parameter  control  law). 
The  fidelity  of  an  in-flight  simulation  is  highly  dependent  on  the 
ability  of  the  variable  stability  system  to  either  complement  the 
natural  stability  derivatives  of  the  host  aircraft  until  they  match 
those  of  the  simulated  vehicle  (response  feedback  mode),  or  to  establish 
a  unity  transfer  function  between  the  responses  of  a  computer  model  of 
the  simulated  vehicle  and  the  responses  of  the  host  aircraft 
(model-following  mode)  (39).  In  both  cases  it  is  necessary  to  have 
detailed  knowledge  of  the  host  aircraft  stability  derivatives  to  perform 
a  satisfactory  simulation  (39,42).  Although  gain  scheduling  is  often 
used  to  alleviate  plant  parameter  variation  problems,  this  can  at  best 
be  considered  open-loop  compensation  which  may  require  extensive  flight 
testing  for  proper  fine  tuning. 


Features 
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Adaptive  control  is  a  promising  approach  to  deal  with  the 
problem  of  maintaining  a  specified  level  of  tracking  performance,  and 
therefore  simulation  fidelity,  throughout  the  operating  envelope.  The 
primary  reason  for  considering  the  use  of  an  adaptive  controller  for  the 
VISTA/F-16  is  the  wide  range  of  dynamic  characteristics  assumed  as  the 
operating  point  changes.  Rather  than  designing  off-line,  fixed  point 
control  strategies,  i.e.  control  strategies  designed  for  a  specific 
operating  condition,  and  attempting  to  chose  the  appropriate  strategy  to 
meet  the  current  conditions,  the  basic  idea  is  to  perform  on-line  system 
identification  and  control  design  simultaneously.  As  the  operating 
conditions  change,  so  does  the  system  being  identified,  and  appropriate 
changes  are  made  to  the  control  law.  Thus,  the  control  system  is 
adapting  to  the  changing  environment.  The  identification  subsystem  uses 
the  inputs  and  outputs  of  the  plant  in  order  to  identify  an  equivalent 
input-output  model  for  the  aircraft  dynamics.  The  identified  parameters 
are  then  used  to  calculate  a  set  of  time-varying  controller  gains,  which 
in  turn  are  used  to  compute  the  current  control  inputs  to  the  plant. 
The  structure  of  such  a  system  is  shown  in  Figure  1-5. 

Much  effort  has  been  devoted  recently  to  the  investigation  of 
adaptive  control  laws  based  on  the  application  of  recursive  parameter 
identification  algorithms.  Positive  results  have  been  reported  in  the 
literature  on  the  use  of  such  adaptive  control  schemes  for  aerospace 
applications  (20,41,48).  New  theoretical  developments  in  the  area  of 
parameter  identification  (19),  and  advances  in  microprocessor  technology 
(1,12,26,47),  make  the  alternative  of  parameter-adaptive  control 


Figure  1-5.  Parameter-adaptive  control  system. 
model  motions;  V  noise  signals. 


increasingly  viable. 

1.4  Problem  Statement 

The  purpose  of  this  thesis  is  to  test  the  effectiveness  of 
digital  parameter-adaptive  control  laws  in  maintaining  tracking 
performance  for  an  in-flight  simulator,  despite  plant  parameter  changes. 
The  use  of  parameter-adaptive  control  laws  will  alleviate  the 
requirement  for  accurate  knowledge  of  the  stability  derivatives  of  the 
host  aircraft  and  offer  the  potential  of  increased  fidelity  for 
in-flight  simulations.  The  design  techniques  of  Professor  Brian  Porter 
(31-38),  are  used  to  develop  longitudinal  control  laws  for  a 
model-following  application.  These  new  techniques  use  on-line. 


recursive,  step-response  matrix  identifiers  to  update  the  control  law 
gains  as  needed  to  account  for  plant  parameter  variations.  A  parameter 
identification  algorithm  recently  developed  by  T.  Hagglund  (19)  is 


chosen  for  this  implementation. 


This  algorithm  allows  for 


identification  during  general  aircraft  maneuvering  by  updating  the 
parameter  estimates  only  when  adequate  input-output  activity  is  present. 

1.5  Scope 

Because  of  budgetary  constraints,  the  initial  VISTA  will  be 
developed  without  side-force  surfaces  and  wingtip  speed  brakes,  in  a 
configuration  such  as  the  one  depicted  in  Figure  1-6.  Since  aerodynamic 
data  is  not  yet  available  for  the  VISTA/F-16,  and  because  of  the 
similarity  between  its  proposed  physical  configuration  and  that  of  the 
Advanced  Technology  Integration  F-16  (AFTI/F-16,  Figure  1-7),  this 
thesis  adresses  the  development  of  adaptive  multivariable  tracker 
control  laws  using  the  AFTI/F-16  as  the  host  aircraft.  The 

investigation  is  limited  to  the  linearized,  longitudinal,  rigid  body 
dynamics  of  the  AFTI/F-16  using  perturbation  equations  of  motion  at  a 
nominal  flight  condition  of  Mach  0.9  at  10,000  ft  MSL. 


This  thesis  accomplishes  the  following  objectives: 

1.  Successful  control  of  the  linear  aircraft  model  including 
actuator  dynamics  with  position  and  rate  limits. 

2.  Determines  the  feasibility  of  using  fast-sampling 
multivariable  tracker  control  laws  in  a  model-following 
configuration  for  in-flight  simulation. 

3.  Successful  implementation  of  state-of-the-art  recursive 
identification  algorithm  for  on-line  tunning  of  control  law 


4.  Assesses  the  improvement  in  tracking  performance  and 
simulation  fidelity  by  using  an  adaptive  system. 

5.  Studies  closed-loop  system  performance  under  simulated 
output  measurement  noise  conditions 

1.6  Overview 

The  material  in  this  thesis  consists  of  a  brief  summary  of 
current  model-following  techniques  presented  in  Chapter  II.  This 
chapter  also  provides  a  description  of  the  multivariable  design  theory 
developed  by  Professor  Brian  Porter  of  the  University  of  Salford, 
England,  presented  here  as  an  alternate  approach  to  the  model-following 
problem.  Chapter  III  then  discusses  the  recursive  technique  used  to 
identify  the  step-response  matrix  elements  required  by  the  control  law. 
Chapter  IV  presents  details  of  the  design  procedure  used  in  the 
development  of  the  control  law  and  recursive  identification  algorithm, 
as  well  as  a  description  of  the  simulation  setup  in  which  the  parameter 
identifier  is  coupled  with  the  control  law  to  form  an  adaptive  system. 
The  capability  of  the  model-following  system  to  adjust  itself  to 
changing  conditions  is  analyzed  through  the  use  of  the  MATRIX^  software 
package  (23).  Chapter  V  presents  simulation  results  and  compares  the 
performance  of  the  adaptive  controller  to  that  of  a  fixed  gain  system. 
Finally,  a  summary  of  the  simulation  results,  as  well  as  conclusions  and 
recommendations  for  further  studies  are  presented  in  Chapter  VI.  Four 
appendices  are  included  as  supplementary  material  to  the  body  of  the 
thesis.  Appendix  A  provides  data  for  the  state  space  and  difference 
equation  models  representing  the  equations  of  motion  of  the  AFTI/F-16 
aircraft.  Appendix  B  details  the  equivalence  between  the  step-response 


matrix  and  the  matrix  coefficient  of  the  vector  difference  equation 
describing  the  input-output  dynamics  of  the  plant.  Appendix  C  covers 
implementation  details  of  the  recursive  parameter  identification 
algorithm,  and  Appendix  D  presents  the  results  of  a  root-locus  analysis 
and  some  of  the  time  responses  used  for  asses ing  the  implications  of 
functional  uncontrollability  on  the  closed-loop  system. 


2.1  Introduction  (39:21-22) 

The  concept  of  using  a  stability  augmentation  system  to  modify 
the  dynamic  characteristics  of  an  airplane  is  not  a  new  one.  The  same 
concept  can  be  used  to  adjust  the  terms  in  the  equations  of  motion  of 
the  variable  stability  airplane  to  match  the  corresponding  terms  in  the 
equations  of  motion  of  the  simulated  airplane.  This  technique  is  known 
as  the  "response  feedback"  approach  to  in-flight  simulation.  A  response 
feedback  variable  stability  system  can  be  described  as  a  generalized 
stability  augmentation  system  which  has  wide  ranges  of  adjustment  so 
that  large  variations  in  airplane  response  characteristics  can  be 
produced  (Figure  2-1).  A  response  feedback  system  operates  by  adding  to 
or  substracting  from  the  airplane's  natural  stability  and  control 
characteristics.  Thus  it  is  necessary  to  know  accurately  the  stability 
and  control  characteristics  of  the  host  airplane  at  all  flight 
conditions  of  interest.  Also,  it  is  difficult  to  calculate  exactly  in 
advance  the  variable  stability  system  gain  settings  which  will 
simultaneously  produce  correct  values  of  the  many  parameters  which 
define  the  aircraft  dynamics.  Thus  in-flight  calibration  of  the 
configurations  that  are  to  be  evaluated  is  generaly  necessary  for 
fine-tuning  of  the  system  in  order  to  match  the  motions  of  the  simulated 
aircraft. 

A  different  approach  to  in-flight  simulation  uses  the  idea  of 
"model-following".  In  this  type  of  system,  the  signals  coming  from  the 
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evaluation  pilot's  cockpit  are  fed  as  inputs  to  a  computer  which  has  in 
it  the  equations  of  motion  of  the  airplane  to  be  simulated.  The  output 
of  this  computer  is  a  set  of  time  histories  describing  the  motion  of  the 
simulated  airplane  to  the  inputs  applied  by  the  pilot.  The  task  of  the 
variable  stability  system  then  is  to  automaticaly  operate  the  control 
surfaces  of  the  host  airplane  in  such  a  way  that  its  motions  (at  the 
pilot's  station)  duplicate  those  of  the  modeled  airplane.  The 

model-following  concept  is  illustrated  in  Figure  2-2.  The 
model-follwing  approach  permits  the  computer  that  defines  the  simulated 
aircraft  to  be  set  up  and  checked  out  on  the  ground  prior  to  flight. 
Although  less  sensitive  to  host  airplane  parameter  variations  than  the 
response  feedback  technique,  the  model  following  approach  still  requires 
an  accurate  knowlwdge  of  the  host  airplane's  stability  derivatives  for 
good  model -fol lowing  performance  as  is  shown  next. 

2.2  Summary  of  current  model -fol lowing  techniques  (39:21-28) 

Under  the  assumptions  presented  in  Etkin  (15:121-189)  the  host 
vehicle  dynamics  at  a  particular  flight  condition  can  be  represented  by 
a  linear  matrix  differential  equation  of  the  form 

'^p(t)  =  Ap  Xp(t)  +  Bp  Up(t)  (2“1) 

and  the  dynamics  of  the  simulated  vehicle  by  the  linear  matrix 
differential  equation 

~  ^m  ^m  ^m(^)  (2-2) 

where  Ap,  A^^,  Bp,  are  matrices  of  stability  derivatives,  /p  and 
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are  the  states,  and  are  control  surface  deflections  and  the 
subscripts  p  and  m  denote  the  plant  (host)  and  model  airplane 
respectively.  Assuming  that  the  simulation  is  started  with  the  host  and 
the  simulated  aircraft  at  the  same  initial  conditions,  a  control  law  for 
exact  model-following  can  then  be  obtained  by  simply  substituting  the 
state,  and  rate  of  change  of  the  state,  of  the  model  into  Eqn  (2-1)  and 
solving  for  the  control  input  i/p(t): 


■^m(^)  ~  '^m(t)  +  Bp  Up(t) 

^^p(t)  =  Bpl  ”  ^p^  ^p  ''^m(^) 

Up{t)  =  Ki  Xm(t)  -  K2  4(t) 


(2-3) 

(2-4) 

(2-5) 


Equations  (2-4)  and  (2-5)  define  the  control  inputs  to  the  host  aircraft 
in  terms  of  the  model  states  and  rate  of  change  of  those  states.  Both 
of  them  are  available  from  the  equations  of  motion  of  the  simulated 
vehicle  that  are  contained  in  the  variable  stability  system's  simulation 
computer  (Figure  2-3). 


PILOT 


'^m(t) 


INPUTS  ;  MODEL  - , 

--  )  AIRCRAFT  X^it) 

:  DYNAMICS  - 


Xn(t) 


Figure  2-3.  Explicit  Model  Following  System 


Equation  (2-4)  also  demonstrates  the  requirement  for  accurate  knowledge 
of  the  stability  derivatives  of  the  host  aircraft. 

To  reduce  the  sensitivity  of  this  control  scheme  to  plant  parameter 
variations,  a  feedback  loop  is  introduced  around  the  plant  with  a  gain 
matrix  Kp  (Figure  2-4). 


Figure  2-4.  Improved  Explicit  Model  Following  System 


The  feedforward  gains  are  now  determined  as  follows: 

Ki  =  Bpl  (2-6) 

K2  =  Bpl  Ap  -  Kp  (2-7) 

The  feedback  around  the  plant  is  arbitrary,  subject  to  the  requirement 
that  the  regulator  loop  must  be  stable.  It  can  be  shown  that  by 
selecting  K2  as  in  Eqn  (2-7),  the  dynamics  of  the  tracking  error 
become 

e(t)  =  (Ap  -  Bp  Kp)  e(t)  (2-8a) 

=  Ag  e(t)  (2-8b) 
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If  Ag  is  stable,  the  error  is  driven  to  zero.  By  selecting  large 
values  for  the  gains  in  the  matrix  Kp  the  regulator  loop  becomes 
increasingly  tight,  making  the  system  less  dependent  on  the  gain  matrix 
and  reducing  its  sensitivity  to  plant  parameter  variations.  This 
reduces  the  requirement  for  accurate  knowledge  of  the  plant's  stability 
derivatives.  In  practice,  concerns  for  closed-loop  stability,  sensor 
noise,  structural  limitations,  etc.  determine  how  high  the  feedback 
gains  can  be,  thus  restricting  the  amount  of  desensitization  (to  plant 
parameter  variations)  that  can  be  provided. 

The  feedback  gains  can  be  selected  based  on  the  maximum  gain 
available  that  will  not  significantly  compromise  closed-loop  stability 
or  result  in  undesirable  control  surface  rates  of  motion  which  could 
place  excessive  demands  on  the  hydraulic  system  used  to  actuate  the 
control  surfaces  of  the  plant.  Usually  they  are  designed  as  constants 
to  reduce  system  complexity,  however,  it  important  to  note  that  the 
effect  on  the  output  of  the  control  surfaces  is  a  function  of  dynamic 
pressure  since  the  actual  control  effectiveness  of  the  surface  is  also  a 
function  of  dynamic  pressure.  Thus  the  aerodynamic  gain  will  vary  with 
flight  condition.  This  increases  the  actual  feedback  gain  as  dynamic 
pressure  increases  and  could  result  in  closed-loop  instabilities.  This 
condition  can  be  minimized  if  the  feedback  gains  are  selected  at  the 
highest  dynamic  pressure  to  be  encountered  in  the  simulation,  although 
this  may  compromise  performance  at  flight  conditions  characterized  by 
low  dynamic  pressure.  An  alternate  solution  is  simply  to  use  gain 
scheduling,  but  this  approach  may  require  extensive  flight  testing  to 
determine  an  appropriate  schedule  for  the  large  number  of  possible 


flight  conditions  likely  to  be  encountered  during  simulations. 

2.3  Porter's  Method 

2.3.1  Fixed  Gain  Control  Law.  One  of  the  many  alternatives  for 
achieving  model  following  is  the  use  of  fast-sampling  multivarible 
control  laws  in  a  configuration  similar  to  a  command  generator/tracker. 
O'"  crol  laws  such  as  the  ones  developed  by  Prof.  B.  Porter  can  be  used 
to  make  the  plant  follow  the  time  histories  generated  by  a  computer 
model  of  the  simulated  vehicle. 

The  basic  longitudinal  equations  of  motion  of  the  host  aircraft  in 
this  thesis  are  assummed  to  be  completely  controllable  and  observable, 
and  described  by  state  and  output  equations  of  the  form 

X{t)  =  A  X(t)  +  B  u{t)  (2-1) 

y(t)  =  C  X{t)  (2-9) 

where  the  subscript  p  has  been  droped  for  rotational  convenience,  and 
A  is  (n  X  n),  B  is  (n  x  m)  and  has  rank  "m",  C  is  (m  x  n),  y(t)  and  u{t) 
are  (m  x  1). 

The  A,  B,  and  C  matrices  are  partitioned  according  to  the  control 
input  matrix  (B)  to  yield  the  following  equations: 


where  x,(t)  is  (p  x  1),  x^{t)  is  (n-p  x  1)  and  B,  is  (n-p  x  m)  with 
rank  "m".  The  equations  defining  XiCt)  are  kinematic  relationships  with 
no  forcing  function  present.  For  the  aircraft  longitudinal  equations  of 
motion  the  kinematic  equation  is  0  =  q,  where  B  is  the  pitch  angle  and 
q  is  pitch  rate.  The  resulting  state  vector  has  the  form: 

^,  =  [  fl  ]  deg  (2-12) 

Xj  =  [  u  a  q  fps,  deg,  deg/sec  (2-13) 

The  plant  is  considered  "regular"  or  "irregular"  depending  on 
whether  or  not  the  first  Markov  parameter  ,  CB,  has  full  rank  equal  to 
"m".  For  "regular"  plants  (rank(CB)  =  m)  with  stable  transmision 
zeros,  the  control  law  is  a  discrete  proportional-plus-integral  (PI) 
output  feedback  control  law  expressed  as 


and 


where 


u(kT)  =  (1/T)  [  K,  e(kT)  +  K,  Z(kT)  ]  (2-14) 

u{t)  =  u(kT)  for  te[kT,  (k+l)T)  (2-15) 


k  =  integer 

T  =  1/f  is  the  sampling  period, 

K,  and  are  (  m  x  m  )  controller  gain  matrices, 

e(kT)  =  r(kT)  -  y(kT)  is  the  (  m  x  1)  error  vector 

'"(kT)  =  /model (I^T)  is  the  reference  vector 

Z(kT)  is  the  digital  integral  of  the  error  vector,  defined  as 


Z[(k+1)T]  =  Z(kT)  +  T  e(kT) 


(2-16) 


For  "irregular"  plants,  a  proportional  plus  integral  plus  derivative 
(PID)  controller  structure  can  be  used  as  described  in  References  (31), 
(32)  and  (33).  The  input-output  configuration  in  this  thesis  provides 
for  a  "regular"  design. 

The  plant  state  and  output  equations  may  be  discretized  for  the 
sampling  period  T  using  the  relationships 
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=  exp 
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(2-17) 


fT  fr  A,,  A,, 

J  A,.  A,,  t 


(2-18) 


The  resulting  sampled  data  state  and  output  equations  for  the  plant  are 


■^,{(1<+1)T}|  rt, .  d),  1  fx,(kT)j  'P 

^j{(k+l)T}  ~  <{>,,  d>,,  x^(kT)  f 


(2-19) 


y(kT)  =  [  C,  C,  ] 


x,{kjy 

x,(kT) 


(2-20) 


The  augmented  closed-loop  state  and  output  equations  for  the  control 
law  of  Eqn  (2-12)  are  given  by  (36) 
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(2-21) 


y(kT)  =  [  0  C,  C,  ] 


Z(kT) 

^,(kT) 

^,(kT) 


(2-22) 


where  f  is  the  sampling  frequency.  A  transformation  is  introduced 
(27)  which  block  diagonalizes  the  state  and  output  equations  so  that 
they  assume  the  form 


x,{(k+l)T}  1 

•A,  0  • 

xa(k+l)T} 

0  A, 

y(kT)  =  [  C,  C, 


where 

x,(kT)  = 

x,(kT)  = 

C,  = 
C.  = 


x,(kT) 

x,(kT) 

B,  • 

r(kT)  (2-23) 
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(2-30) 
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(2-31) 

1 

(2-32) 

As  the  sampling  frequency  is  increased  (f  -»  »),  the  closed-loop 
transfer  function  assumes  the  asymptotic  form  (37) 

r(2)  =  r,(2)  +  r,(2)  (2-33) 

where  2  is  the  discrete  transform  operator  and 

r,(2)  =  C,(zIn-ln-TA„)-lTB„  (2-34) 

r,(2)  =  C,(zIn,-Im-Aj-lB,  (2-35) 

with 


Ao  = 


K7IK,  0 

A,.C7lK7lK,  A.,-A,,C7lC, 


(2-36) 


Bo  = 


A,,C7l 


(2-37) 


A,  =  -B,K,C, 


(2-38) 


r,(2)  and  rj(2)  are  the  slow  and  fast  mode  transfer  functions 
respectively.  The  slow  modes  can  be  grouped  into  two  sets  Z,  and  Z. 
and  are  given  by 
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Z,  =  {zeC:  TK.'Ik.i  =  0} 


(2-39) 


and 

Z,  =  {zeC:  izIn.n,-In.ni-TA,,+  =  0}  (2-40) 

The  fast  modes  are  given  by 

Z_,  =  {zeC:  izln,-lm+  C,B/.  i  =  0}  (2-41) 

Because  of  the  form  of  ,  B„  ,  and  "C,  ,  the  eigenvalues  of  A„ 
are  uncontrollable  or  unobservable.  Thus,  as  the  sampling  frequency 
increases,  the  slow  transfer  function  asymptotically  approaches  zp’^n  and 
the  overall  system  transfer  function  contains  only  the  tast  modes,  as 
given  by  rj(z)  which  can  be  put  in  the  equivalent  form 

r(z)  =  r,(z)  =  (zlni-lni+  C,B,K,  )-lC,B,K,  (2-42) 

The  controller  matrices  K,  and  are  then  given  by 

K,=  [C,BJ-1e  (2-43) 

K,=  pK,  (2-44) 

where  p  is  any  positive  scalar  greater  than  zero,  and  E  is  a 

diagonal  tuning  matrix.  Both  p  and  E  are  chosen  by  the  designer  to 
achieve  the  desired  tracking  characteristics. 

Although  the  method  just  outlined  provides  robust  control 
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characteristics,  some  degradation  of  performance  is  inevitable  when 
F  ced  with  large  parameter  variations  in  the  plant.  To  avoid  this 
performance  degradation,  adaptive  control  techniques  may  be  used  to 
adjust  the  control  law  parameters  as  necessary  and  maintain  tracking 
performance. 


2.3.2  Adaptive  Control  Law.  An  alternate  way  of  defining  the 
control  law  gain  matrices  of  Eqns  (2-43)  and  (2-44)  is  by  the  use  of  the 
step-response  matrix  (34)  defined  as 


H(T)  =  J  C  exp(A(t-T))  B  d7 


(2-45) 


For  small  sampling  periods  H(T)  s  TCB  ,  and  the  control  law  can  be 


expressed  as 


where 


//(kT)  =  K,e(kT)  +  K,Z(kT) 


K,  =  h-1(T)i; 
K,  =  pK, 


(2-46) 


(2-47) 

(2-48) 


The  significance  of  using  H(T)  is  that  it  can  be  obtained  from 
real-time  input-output  measurements  to  reflect  the  current 
characteristics  of  the  plant  (31,32),  thus  forming  the  basis  for  an 
adaptive  system.  In  order  to  do  this,  an  autoregressive  difference 
equation  is  used  to  represent  the  open-loop  dynamics  of  the  host 
aircraft.  This  is  accomplished  in  the  next  chapter. 
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III.  Recursive  Identification 


3.1  Introduction 

In  recent  years  numerous  new  developments  for 
parameter-adaptive  control  algorithms  based  on  parameter  estimation 
techniques  have  taken  place.  Recursive  Least  Squares  identification 
algorithms  play  a  crucial  role  for  many  problems,  not  only  in  adaptive 
control,  but  also  for  adaptive  signal  processing,  and  for  general  model 
building  and  monitoring  problems.  Today,  with  the  aid  of  an  on-line 
digital  computer,  it  is  possible  to  use  an  empirical  model  to  represent 
a  non-linear  plant  in  a  sequential  manner.  Specifically,  a  low  order 
linear  model  can  be  selected  to  represent  the  dynamics  of  the  system. 
By  use  of  the  proper  estimation  algorithm,  the  computer  is  able  to 
provide  updated  parameter  estimates  for  the  empirical  model  at  each 
sample  instant.  This  quasi-linearization  procedure  permits  the  use  of 
advanced  control  techniques  such  as  the  self-adaptive  strategy,  which 
uses  the  information  about  the  empirical  model  to  update  control  law 
parameters,  thus  allowing  improved  performance  over  a  wide  range  of 
system  parameter  variations.  This  chapter  presents  a  short  description 
of  the  procedure  and  algorithms  used  to  carry  out  the  identification 
part  of  the  adaptation  mechanism. 

The  identification  problem  considered  in  this  thesis  is  the 
following: 

Given  a  discrete  input-output  time  history  of  the  airplane 
longitudinal  dynamics,  determine  the  parameters  of  a  suitable  linear 


model  that  best  describes  the  characteristics  of  the  airplane  (according 
to  a  given  criterion)  at  its  current  operating  condition.  Furthermore, 
the  identification  must  be  done  "on-line"  with  a  minimal  amount  of  data 
storage  (recursive  algorithm  implementation),  and  using  only  the  normal 
control  inputs  to  the  plant. 

Section  3.2  of  this  chapter  describes  the  process  of  selecting  the 
structure  of  the  input-output  model  used  in  the  identification 
procedure.  Section  3.3  then  sumarizes  the  commonly  used  Recursive  Least 
Squares  identification  algorithm,  its  characteristics,  and  some  of  its 
limitations.  This  leads  the  way  to  Section  3.4  which  sumarizes  a 
modified  version  of  the  Least  Squres  algorithm  developed  by  Hagglund 
(19)  which  overcomes  some  of  the  limitations  of  the  basic  Least  Squares 
approach.  For  this  investigation,  it  will  be  assumed  that  the 
computational  time  involved  in  the  execution  of  the  control  and 
identification  algorithms  is  sufficiently  short  (compared  to  the 
sampling  time,  and  variation  of  plant  parameters)  as  to  introduce 
negligible  time  delay  effects.  This  assumption  is  justified  by  the 
availability  of  sophisticated  microproccesors  such  as  the  ones  described 
in  References  (12),  (26),  and  (47)  which  may  allow  for  fast  execution  of 
these  algorithms. 

3.2  Input-output  Model  Specification 

The  first  step  of  the  identification  problem  is  the 
specification  of  an  algebraic  structure  between  the  input  and  output 
variables  to  describe  the  behavior  of  the  plant.  This  model  is 


postulated  a  priori  and  can  be  based  upon  knowledge  of  the  underlying 
theory  governing  the  dynamics  of  the  plant.  The  discrete  nature  of  the 


control  laws  being  considered  here,  in  which  input-output  data  is 
gathered  by  the  on-line  computer  for  the  control  process,  suggests  that 
the  most  convenient  model  is  in  the  form  of  difference  equations. 

One  of  the  posibilities  for  the  derivation  of  a  difference  equation 
model  for  the  open-loop  longitudinal  dynamics  of  the  host  airplane  is 
based  on  the  discrete  state  and  output  relationships  given  in  Eqns 
(2-17)  and  (2-18).  These  equations  can  provide  the  desired  input-output 
relationship  by  simply  taking  their  1  Transform  and  obtaining  a  transfer 
function  model  (22)  as  follows: 

X  {(k+l)T}  =  $  X(kT)  +  ^  u{kl)  (3-1) 

y(kT)  =  C  ;f(kT)  (3-2) 

Y(z)  =  [  C  [zl  -  <t]-l  y  ]  U(z)  =  G(z)  U(z)  (3-3) 

The  m  X  1  vector  Y(z),  the  m  x  m  matrix  G(z),  and  the  m  x  1  vector 
U(z)  are  given  by: 


Y(z) 


Yi(z) 

Y2(z) 

Ym(z) 


(3-4) 
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1  Gmi(z) 

Gm2(z) 

•••  Gfjini(z) 

(3-5) 


where  G^jCz)  is  the  transfer  function  relating  the  output  Y-j  to  the 
control  input  Uj  and  is  of  the  form 

bi  zw  +  b2  +  *••  +  by^  z  +  b^+i 

Gij(z)  = - ^ -  (w  <  n)  (3-7) 

zn  +  ai  zn-1  +  ...  +  an-i  z  +  ap 

By  dividing  each  numerator  and  denominator  in  G(z)  by  z'^  the 
transfer  matrix  G(z)  is  transformed  into  the  delay  operator  form 
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bi  zw-n  +  b2  zW-n-1  +  ...  +  z""^ 

1  +  ai  z-1  +  ...  +  an-i  z'^+l  +  ap  z-n 


(w  <  n)  (3-8) 


At  this  point,  the  input-output  relationship  of  the  open-loop  plant  may 
be  obtained  by  grouping  the  coeficients  with  the  same  amount  of  delay  in 
every  transfer  function  G-jj(z),  and  taking  the  inverse  Z  transform  to 
yield  the  following  N^h  order  vector  difference  equation  model 


y(kT)  =  Bi  o{(k-l)T}  -  Ai  y{(k-l)T}  +  .  .  . 

+  Bn  a{(k-N)T}  -  An  y{(k-N)T}  +  c(kT)  (3-9) 


or  equivalently 


y(kT)  =  TT(kT)  0  +  c(kT) 


e  Rni 


(3-10) 


(here 

t(kT)  is  an  equation  error  term  assumed  to  be  a  zero  mean 
Gaussian  white-noise  vector  with  elements  of  variance 
added  to  account  for  modeling  errors  that  may  arise  from 
applying  a  linear  model  structure  to  a  possibly  nonlinear 
plant,  or  mismatch  between  the  real  system  order  "n"  and  the 
difference  equation  order  "N" 

T^CkT)  e  is  a  matrix  of  past  values  of  {y(kT)}  and  {u(kT)] 

the  matrices  A-j  e  (  i  =  1,2,...,N), 

(  1  =  1,2,...,N)  and  the  vector  o  e  rT><1  are  the  parameters 
of  the  Nth  order  difference  equation 

By  defintion  of  the  step-response  matrix  it  can  be  shown  (Appendix  B) 
that 

H(T)  =  Bi  s  TCB  (3-11) 


thus,  by  identifying  the  matrix  coefficient  Bj  in  real-time  from 
input-output  data,  and  invoking  the  certainty  equivalence  principle  (6), 
updated  step-response  matrix  estimates  can  be  provided  for  the  control 
law  design  calculations  (Eqn  (2-47))  (31,32,34). 

It  is  important  to  note  that  all  the  parameters  of  Eqn  (3-9)  must  be 
identified  to  obtain  the  step-response  matrix  denoted  by  Ideally, 
in  the  case  a  transfer  matrix  model  reduced  to  the  same  denominator,  the 


order  "N"  should  be  the  same  as  the  discretized  plant  model  order  "n" 
(Eqns  (2-17)  and  (2-18))  to  avoid  biases  in  the  estimates  resulting  from 
trying  to  fit  a  reduced  order  transfer  function  model  to  the  plant 
dynamics  (  this  assumes  of  course  that  a  model  of  order  "n"  is  true  ). 

A  model  in  which  "N"  is  large  however  may  represent  a  problem,  since  the 
computational  effort  and  the  convergence  time  of  the  identification 
algorithm  depends,  among  other  factors,  on  the  number  of  parameters  to 
be  identified.  It  is  important  therefore  to  keep  the  number  of 
identified  parameters  to  a  minimun.  This  fact  seems  to  oppose  the 
desire  of  having  "N"  equal  to  "n",  favoring  the  use  of  reduced  order 
models.  A  reduced  order  model  may  be  appropriate  for  certain  situations 
depending  upon  the  tracking  performance  required  and  desired  gain 
margins.  Compensation  for  biases  in  the  step-response  matrix  can  be 
introduced  by  increasing  the  overall  loop  gain  of  the  system.  However, 
it  is  highly  desirable  to  achieve  the  desired  tracking  behavior  with  the 
lowest  amount  of  gain  possible  in  order  to  provide  satisfactory  gain 
margins  in  the  system.  A  solution  is  then  needed  that  pe’^mits  the 
reduction  of  biases  in  the  parameter  estimates  associated  with  reduced 
order  modeling,  while  also  providing  for  a  small  amount  of  parameters  to 
be  identified.  It  is  recognized  that  a  diference  equation  model  derived 
from  a  transfer  matrix  does  not  provide  for  a  minimal  number  of 
parameters  to  represent  a  system  (17,18),  however,  it  facilitates  the 
identification  procedure  by  providing  for  decoupled  equation  error  terms 
which  in  turn  allow  for  a  simple  scalar  measurement  update  procedure 
(10,  30,  43,  Appendix  C). 

Since  only  (=  H(T))  is  required  for  the  control  law  calculation 
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of  Eqn  (2-47),  a  possible  approach  to  solve  the  aforementioned  problems 
is  to  identify  only  the  elements  of  leaving  the  rest  of  the 
parameters  in  the  full  order  model  (  N  =  n  )  fixed.  It  is  assummed  for 
this  investigation  that  the  appropriate  parameters  for  the  fixed  portion 
of  the  model  of  Eqns  (3-9)  and  (3-10)  (representative  of  the  current 
flight  condition)  are  available  for  use  in  calculating  the  plant's 
output  predictions  and  residuals  needed  in  the  identification  algorithm. 
This  approach  is  justified  in  light  of  the  results  obtained  in 
references  (3),  (21),  and  (46),  where  multiple  models  are  employed  for 
identification  based  on  a-priori  data.  This  a-priori  data  can  be 
provided  from  information  gathered  by  stability  and  control  derivative 
prediction  methods,  wind  tunnel  tests,  and  preliminary  flight  test  data. 
Selection  of  the  best  performing  model  could  then  be  based  upon  on-line 
calculation  of  a  probability  of  correctness  associated  with  each  model. 
Furthermore,  the  same  input-output  data  used  in  updating  the  estimates 
of  the  step  response-matrix  elements  can  be  used  to  update  a  full  order 
model  that  may  be  evaluated  on  a  background  proccesing  mode.  This 
information  can  eventually  be  used  in  refining  the  models  selected 
a-priori  for  improved  performance.  The  computational  effort  involved  in 
the  multiple  model  approach  is  easily  accomodated  by  pararallel 
processing  techniques  and  should  not  produce  any  significant 
computational  time  delays. 

The  procedure  of  identifying  only  causes  the  following 

partition  of  Eqn  (3-10) 

y(kT)  =  TT(kT)  o'  +  Q(kT)  +  v(kT)  t  R~i  (3-12) 

in  which  the  known  or  fixed  parameters  in  o  have  teen  deleted  to  form 


e'.  This  is  also  done  with  the  corresponding  elements  in  T(kT).  These 
non  identified  parameters  and  appropriate  delayed  measurements  are  then 
collected  into  the  vector  o(kT).  The  vector  e'  thus  contains  only 
the  elements  of  and  'f'(kT)  the  corresponding  elements  of  T(kT). 
Additional  details  are  given  in  Appendix  C. 

Having  settled  upon  a  model  structure  to  represent  the  plant 
dynamics,  it  is  now  necessary  to  choose  a  recursive  algorithm  to 
identify  the  required  elements  of  the  step-response  matrix.  Sections 
3.3  and  3.4  are  devoted  to  this  subject. 


Given  a  linear  model  of  the  structure  described  by  Eqn  (3-9), 
our  goal  is  to  find  an  estimate  of  the  parameter  vector  0  such  that  it 
minimizes  the  following  weighted  quadratic  cost  criterion  on  the 
parameter  vector  estimate 

N 

0=1  {  [  y(kT)  -  TT(kT)  0  ]2  }  (3-13) 
k=0 

where  N  represents  the  total  number  of  measurements,  and  the  quantity 
within  the  brackets  is  the  error  between  the  measurement  of  the  actual 
output  of  the  system  and  the  predicted  output  based  on  the  current 
estimate  of  0.  Differentiating  the  cost  function  with  respect  to  the 

1  This  discussion  pertains  to  a  single-input  single-output 
difference  equation  model  although  the  multiple-input  multiple-output 
case  can  be  easily  accomodated  by  performing  a  scalar  measurement  update 
procedure  [  (10),  (28),  Appendix  C  ]. 


parameter  vector  (0)  and  after  some  derivations  (28),  the  recursive 
least-squares  (RLS)  identification  algorithm  is  obtained 


RCkT)  =  R{(k-1)T}  +  -4t-  [  'r(kT)  T(kT)T  -  R{(k-1)T}  ]  (3-14) 


Y(kT)  =  y(kT)  -  T(kT)T  G{(k-1)T} 


0(kT)  =  0{(k-l)T}  +  — R(kT)-l  T(kT)  T(kT) 


(3-15) 


(3-16) 


where  R(kT)  is  the  information  matrix,  and  T(kT)  is  the  prediction 
error.  Though  these  equations  are  recursive,  they  are  not  well  suited 
for  on-line  implementation  since  at  each  step  a  matrix  inverse  is 
required.  However,  taking  the  advantage  of  the  fact  that  the  update  to 
the  information  matrix  (  R(kT)  )  is  of  rank  one,  the  matrix  inversion 
lemma  can  be  used  to  exchange  the  matrix  inverse  for  a  scalar  division 
in  the  propagation  of  the  parameter  covariance  matrix  as  shown  in  Eqn 
(3-17) 


P(kT)  =  P{(k-1)T}  - 


1  +  TT(kT)  P{(k-1)T}  T(kT) 


(3-17) 


where 


P(kT)  =  —  R"l(kT)  is  the  parameter  covariance  matrix 
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The  parameter  update  equation  is 


Q(kT)  =  0{(k-l)T}  +  P(kT)  T(kT)  T(kT) 

=  G{(k-1)T}  -f  P{(k-1)T}  T(kT)  T(kT) 

1  +  TT(kT)  P{{k-1)T}  T(kT) 


(3-18a) 

(3-18b) 


The  previous  equations  weight  each  measurement  equally  and  assume  that 
the  parameters  and  measurement  noise  levels  are  constant  over  time, 
though  the  estimate  o(kT)  is  certainly  a  time-varying  quantity.  The 
interpretation  given  to  G(kT)  is  that  of  the  best  estimate  of  the 
constant  parameter  vector  given  all  the  past  information  up  to  time  kT. 

A  more  realistic  approach  allows  for  time  varying  noise  levels 
and/or  parameter  vector.  Such  is  the  case  in  the  weighted  least-squares 
approach  where  the  cost  function  is  modified  as  follows: 


k — 0 

where  the  weights  W(kT)  are  to  be  selected  in  such  a  way  as  to 
indicate  the  degree  of  confidence  that  can  be  placed  on  the  individual 
measurements,  or  equivalently,  the  amount  of  uncertainty  associated  with 
a  particular  measurement.  A  desirable  choice  for  the  weights  W(kT) 
would  be  the  variances  of  the  corresponding  measurements.  A  key  problem 
in  identification  of  time-varying  systems  is,  however,  the  lack  of 
knowledge  of  these  variances  (19:21). 


Assuming  that  an  estimate  of  the  noise  variance  is  available,  the 
RLS  algorithm  is  given  by: 


0(kT)  =  0{(k-l)T}  +  - p{(k-l)T}  T(kT)  t(kT) _  ^2_20) 

v(kT)  +  TT(kT)  P{(k-1)T}  T(kT) 

P(kT)  =  P{(k-1)T}  -  P{(k-1)T|  T(kT)_TT(kT)  P((k-1)T?, 

v(kT)  +  TT(kT)  P{(k-1)T}  T(kT) 


where  v(kT)  =  W(kT)  is  an  estimate  of  the  noise  variance  at  time  kT. 
Although  this  algorithm  still  assumes  that  the  parameters  are  constant, 
it  allows  for  time  varying  noise  characteristics. 

In  order  to  allow  for  time-varying  parameters,  past  information  must 
somehow  be  deweighted  and  greater  emphasis  placed  on  more  recent 
information. 

A  common  and  simple  way  to  accomplish  this  goal  is  by  modifying  the 
cost  function  to  include  a  "forgetting  factor"  to  cause  an  exponential 
deweighting  of  data  and  place  more  emphasis  on  recent  measurements.  The 
cost  function  then  becomes: 

^  1 

J  ^  0  ]2  xN  -  k  }  (3-22) 
k —0 

where  0  <  x  <  1  is  the  forgetting  factor.  This  cost  function  leads  to 
one  of  the  most  commonly  used  version  of  the  RLS  algorithm: 
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G(kT)  =  e{(k-l)T}  + 


P{(k-1)T}  T(kT' 


X  v(kT)  +  TT(kT)  P{(k-1)T}  T(kT) 


P(kT) 


P{(k-1)T}  - 


X  v(kT)  +  TT(kT)  P{(k-1)T}  T(kT) 


(3-23) 


(3-24) 


The  use  of  the  RLS  algorithm  is  subject  to  a  number  of  preconditions 
and  limitations.  A  principal  assumption  is  that  the  time  variations  of 
the  parameters,  and  noise  level,  are  slow  and/or  seldom  compared  with 
the  time  constants  of  the  system  (19:25);  in  other  words,  that  large 
step-like  changes  in  the  parameters  may  not  occur  frequently. 

An  important  requirement  is  that  of  "presistency  of  excitation"  of 
the  input  data.  To  estimate  a  model  of  the  plant  dynamics  that  is 

satisfactory,  the  input  signals  to  the  plant  must  have  sufficient  energy 
and  rich  frequency  content  within  the  control  bandwidth  of  interest. 
The  conditions  on  persistency  of  excitation  are  related  to  the 
complexity  of  the  estimated  model.  This  implies  that  the  requirements 
on  the  input  signal  become  more  severe  if  the  model  order  (or  the  number 
of  identified  parameters)  is  increased  (45).  Since  the  input  signals  to 
the  plant  are  generated  by  feedback  there  is  no  guarantee  that  these 
signals  will  be  persistently  exciting.  On  the  contrary,  good  regulation 
may  give  a  poor  excitation  (5,  45)  since  the  algorithms  usually  extract 
information  from  the  perturbations  of  the  input-output  signals  around 
the  nominal  set-points.  In  that  case  the  perturbations  can  be  small, 
unless  the  command  signals  themselves  are  dynamic  enough  to  provide  the 
necessary  control  activity.  This  however  is  not  always  possible  or 
practical . 


To  track  parameter  variations  it  is  necessary  to  discount  old  data 
by  means  of  the  forgeting  factor  in  Eqns  (3-22)  and  (3-24).  The  choice 
of  a  suitable  forgetting  factor  is  usually  the  result  of  a  compromise 
between  fast  adaptation  and  high  stationary  accuracy  of  the  parameter 
estimates.  If  the  forgetting  factor  is  small,  old  data  is  discounted 
quickly  and  the  estimated  parameters  will  converge  rapidly  towards  the 
new  values.  However,  the  accuracy  of  the  estimates  will  decrease.  A 
high  value  for  the  forgrtting  factor  (  x  1  )  on  the  other  hand  will 
make  it  impossible  to  track  rapid  parameter  variations  (Figure  3-1). 
The  final  selection  of  a  forgeting  factor  then  reflects  the  compromise 
between  the  demands  on  convergence  rate  and  long  term  quality  of  the 
parameter  estimates.  This  trade-off  is  not  always  satisfactory  (19:17). 
It  is  desirable  to  discount  old  data  quickly  when  the  plant  is  changing 
or  has  just  changed  and  to  discount  data  slowly  when  the  parameters  are 
relatively  constant. 

Another  problem  that  may  occur  when  a  constant  forgetting  factor 
less  than  one  is  used  is  that  of  instabilities  in  the  estimation 
algorithm  commonly  refered  to  as  "estimator  wind-up"  (2).  Exponential 
forgetting  works  well  only  if  the  plant  is  properly  excited  all  the 
time.  This  may  not  be  the  case  if  the  main  source  of  excitation  are 
changes  in  the  set-points.  Then  there  may  be  long  periods  with  no 
excitation,  the  estimator  will  continue  to  discount  old  information  and 
the  uncertainty  of  the  parameters  will  grow.  The  effect  can  be  seen 
analytically.  It  can  be  shown  that  Eqn  (3-24)  is  equivalent  to: 

P(kT)  =  [  v-l(kT)  T(kT)  T‘''(kT)  +  x  P-l{(k-l)T}  ]-l  (3-25) 
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If  the  plant  is  poorly  excited,  the  input  u(kT)  and  the  output  y(l<T) 
are  small.  Since  the  components  of  the  vector  r(kT)  are  delayed 
values  of  the  input  and  output  the  vector  T(kT)  will  be  also  small. 
In  the  extreme  case  when  T(kT)  is  zero  Eqn  (3-25)  leads  to 

P(kT)  =  (3_26) 

indicating  that  the  matrixx  P(kT)  grows  exponentially.  If  there  is  no 
excitation  for  a  long  period  of  timje  then  P(kT)  may  become  very  large 
and  the  estimator  becomes  unstable.  Since  P(kT)  is  the  gain  in  Eqns 
(3-18a)  and  (3-23),  small  prediction  errors  can  produce  large  variations 
in  the  parameter  estimates.  This  in  turn  may  drive  the  closed-loop 
system  unstable. 

Another  way  of  looking  at  the  problem  is  from  a  point  of  view  that 
considers  how  the  information  is  processed  in  the  RLS  algorithm. 
Equation  (3-24)  can  also  be  expressed  as: 


P-l(kT)  =  P‘l{(k-1)T}  -  (1-x)  P-l{(k-l)T}  +  v-l(kT)  T(kT)  T(kT)T  (3-27) 

where  P~l(kT)  is  the  information  matrix.  From  Eqn  (3-27)  it  can  be  seen 
that  the  first  term  represents  the  old  information,  the  second  term 
represents  the  information  taken  out,  and  the  third  term  represents  the 
information  that  is  added  by  the  measurement.  Careful  examination  of 
Eqn  (3-27)  shows  that  data  is  discounted  in  all  directions  but  the 
measurement  is  bringing  information  in  only  one  direction  (the  first  two 
terms  of  Eqn  (3-27)  are  full  rank  while  the  third  term  is  of  rank  one). 
Because  of  this,  if  persistency  of  excitation  is  not  present,  more 


information  will  be  discounted  than  brought  in  and  eventually  P'l(kT) 
will  go  to  zero  thereby  causing  instabilities  in  the  algorithm. 

There  are  several  ways  to  avoid  estimator  wind-up.  One  is  obviously 
making  sure  that  the  plant  is  properly  excited  before  performing  the 
identification.  The  condition  for  persistent  excitation  can  be 
monitored  and  perturbations  may  be  introduced  if  the  excitation  is  poor. 
In  cases  where  perturbations  are  not  convenient,  the  identification 
process  can  be  discontinued  temporarily.  Supe-^visory  loops  ((5,  24), 
Figure  3-2)  are  often  employed  to  perform  these  functions.  The  price 
paid  for  this  safeguard  is  the  extra  logic  and  data  storage  required  to 
perform  the  appropriate  checks. 

Another  possibility  of  dealing  with  the  changing  demands  for  data 
discounting  is  the  use  of  time-variable  forgetting  factors  (16)  to 
adjust  the  amount  of  discounting  automatically.  The  use  of  variable 
forgetting  factors,  however,  is  often  heavily  dependent  on  the 
assumptions  imposed  on  parameter  and  noise  level  variations.  Typically, 
the  forgetting  factor  is  adjusted  in  a  manner  inversely  related  to  the 
prediction  error  or  an  estimate  of  its  variance.  The  prediction  error 
can  be  thought  of  as  being  composed  of  two  independent  components  as 
indicated  below: 

>  =  y(kT)  -  ,T(kT)  c(kT)  =  r,(kT)  +  .p,(kT)  (3-28) 

where  model  error  and  the  measurement  noise. 
In  the  least-squares  method,  each  reasunement  is  weighted  according  to 
an  estimate  of  its  uncertainty  (  Eqn  3-1?  ),  which  can  be  expressed  as 


Figure  3-2.  Parameter-adaptive  Control  system  with 
supervision  functions.  (24) 


fol lows: 


c2(kT)  =  c^(kT)  +  c^(kT)  (3-29) 

where  ^^(kT)  is  the  estimate  of  the  prediction  error  variance  at  time 
kT,  c^(kT)  and  c^(kT)  are  the  model  error  variance  and  measurement 
noise  variance  at  time  kT  respectively.  In  the  case  of  small  and 
constant  noise  levels,  assuming  that  the  plant  is  prope’"ly  excited,  any 
changes  in  the  plant  will  be  reflected  in  the  prediction  error  and  the 
forgetting  factor  will  be  adjusted  accordingly  for  imp-'oved  convergence 
to  the  new  parameter  values.  The  problem  with  this  technique  is  that  an 
increase  of  the  noise  level  will  in  most  algorithms  be  interpreted  as  a 
variation  of  the  parameters,  i.e.  an  increase  in  C[^,(kT).  This  action 
is  a  serious  mistake  (19:28)  since  the  algorithm  will  discount  old  (and 
possibly  good)  information  to  favor  new  measurements  with  poor 
information  content. 
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In  order  to  avoid  these  problems,  and  the  required  assumptions  about 
how  the  parameters  and  noise  level  vary,  the  approach  taken  in  this 
thesis  is  to  use  a  modified  RLS  algorithm  developed  by  Hagglund  (19). 
This  algorithm  automatically  discounts  old  data  depending  on  the  amount 
of  incoming  information  and  updates  the  parameter  estimates  only  in  the 
direction  where  new  information  is  entering.  Section  3.4  discusses  this 
algorithm  in  more  detail. 


Haqqlund's  Algorithm  (19) 


To  account  for  time  varying  parameters  and  to  remedy  the 
weighting  problem  discussed  in  the  previous  section,  a  new  principle  of 
forgetting  old  data  was  presented  by  hagglund  (19)  and  used  in  a 
modified  RLS  algorithm.  This  section  summarizes  the  development  of 
Hagglund's  algorithm. 

As  discussed  earlier,  the  purpose  of  the  least-squares  estimator  is 
to  provide  rasonably  accurate  estimates  of  the  parameter  vector  e(kT). 
To  accomplish  this  objective,  Hagglund  approaches  the  information 
weighting  problem  from  a  different  point  of  view.  Typically,  the 
information  about  the  uncertainties  in  the  mieasurements  is  rather  pocr. 
Thus  instead  of  using  assu~cticns  on  how  the  parameters  and  noise  level 
vary,  the  information  weignting  is  handUj  by  relatir'g  the  accuracy  of 
the  parameter  estimates  i.e.  to  the  amount  cf  infcrmation  availade,  anc 
the  incoming  infcrmation.  Hagglum  developed  his  algorithm  adhe-'ing  tc 
the  following  principle: 
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"Discount  old  data  in  such  a  way  that  a  constant 
desired  amount  of  information  is  retained,  if 
the  parameters  are  constant  (19:70)". 

Ttie  amount  of  information  is  signified  by  the  inverse  of  the  P  matrix. 
The  goal  of  the  estimator  developed  by  Hagglund  is  to  weight  the 
incoming  data  so  that  the  covariance  matrix  becomes  proportional  to  the 
identity  matrix.  The  diagonal  elements  of  P(kT)  may  be  interpreted  as 
approximations  of  the  variances  of  the  corresponding  parameters.  The 
weights  W(kT)  in  Eqn  (3-19)  are  chosen  so  that  these  variances  get  a 
desired  value.  As  a  result,  the  amount  of  the  data  included  in 
calculating  an  estimate  of  the  parameter  vector  is  dependent  on  the 
information  the  estimator  is  receiving.  If  there  is  no  information 
coming  in  nothing  will  be  forgotten.  If  the  incoming  information  is 
small  the  convergence  time  will  be  long.  However,  if  the  information 
content  in  the  data  is  large,  old  measurem,ents  will  be  discounted 
quickly  so  that  fast  parameter  adaptation  can  occur. 

If  data  is  to  be  discounted  according  to  the  new  principle,  the 
information  discounted  must  be  the  same  as  the  information  brought  in  by 
the  new  measurements.  Recalling  from  Eqn  (3-27),  the  new  information  is 
proportional  to  T(kT)  T'^(kT).  It  ray  be  said  that  the  new  information 
is  coming  in  she  direction  cf  i'v^T).  Thus  old  information  is  to  be 
discounted  in  the  same  direction.  In  terms  of  the  information  matrix 
this  is  expressed  as: 


p-l(kT)  =  P-^^(k-l)T}  +  v-l(.,T)  -|(kT)r^(kT) 

-  .(kT)  T(kT)YT(kT)  (3-30a) 


p-l(kT)  =  P-l{(k-l)T}  +  [  v-l(kT)  -  o(kT)  ]  T(kT)  T^(kT)  (3-30b) 


where  a(kT)  is  a  discounting  factor.  Note  that  in  the  updating  of  Eqn 
(3-30),  information  is  removed  only  in  the  direction  that  information  is 
added.  This  is  in  contrast  to  the  earlier  discusson  of  the 
least-squares  algorithm  with  exponential  discounting  in  Eqn  (3-27)  where 
information  was  removed  in  all  directions  due  to  the  second  term  being 
full  rank.  The  information  matrix  equation  (3-30)  gives  way  to  the  new 
method  of  updating  the  covariance  matrix  in  Hagglund's  modified  RLS 
algorithm: 


P(kT)  =  P{(k-1)T} 


. PUk-]).TJ  .KkT),  -^(1,1)  F1(I<-1)T1 _ 

[v-l(kT)  -  a(kT)]-l  +  •,J(kT)  P{(k-1)T}  i(kT) 


(3-31) 


Sirce  the  form  of  the  covariance  matrix  update  has  changed,  the 
equation  for  updating  the  parameter  estimates  will  also  change.  Using 
Eqn  (3-31)  along  with  the  basic  definition  of  the  parameter  update 
equaticn  in  the  least-squares  algorithm,  Hagglund  derives  the  follcwing 
param,eter  estimate  update  equation: 


;(kT)  =  c{(k-l)T} 


P((k-1)T}  V(kT)  M/:T) 


+  -  S  .  .  .  .  .  _  .  .  , -  _  _ 

v(kT)  +  YT(k.T)  P{(k-1)T}  ukT)  [l-4kT)v(kT)] 


-Za) 


=  C{(k-1)T}  +  ,  P(kT)  v(kT:  .^kT) 


(•-330) 


It  remains  to  be  shown  the  choice  of  an  appropriate  discounting 
factor  o(kT)  in  Eqns  (3-31)  and  (3-32a).  Equation  (3-30)  shows  that 
a(kT)  must  be  positive  or  information  would  be  added  instead  of 
removed.  Also,  if  o(kT)  is  too  large  the  covariance  matrix  could 
become  non  positive  definite.  These  and  other  considerations  establish 
the  need  for  a  set  bounds  to  which  the  discounting  factor  must  be 
limited.  In  deriving  a  set  of  bounds  on  ft(kT),  Hagglund  performed  a 
stability  investigation  which  included  showing  that  a  proper  choice  of 
(.(kT)  would  ensure  that  the  covariance  matrixx  remained  positive 
definite.  The  stability  investigation  yielded  the  bounds  on  Q(kT)  and 
were  chosen  such  that 

0  <  a(kT)  <  - - -  (3-33) 

TT(kT)  P{(k-1)T}  Y(kT) 

Furthermore,  in  order  to  obtain  a  diagonal  P-matrix  of  the  form  a*  I 
rthere  a  is  the  desired  variance  of  the  parameter  estimates,  Hagglund 
shows  that  n(kT)  must  be  selected  so  that 


T'UT)  Pi(k-1)T|  P(kT)  Pl(k-l)T}_Y(kT) 
YT(kT)  P{(k-1)T]  P{(k-1)T}  Y(kT) 


=  a 


(3-34) 


Substituting  Eqn  (3-31)  into  Eqn  (3-34)  yield  the  following  desired 
value  of  <>(kT) 


.(kT)  =  v-l(kT) 


r.(kT)  Y"(kT)  F;:k-:;T;  -  i 


(3-3F) 


6jj(kT)  = - ^ -  a  (3-36) 

T^(kT)  P2{(k-1)T}  T(kT)  TT(kT)  p2{(k-l)T}  T(kT) 


To  gain  seme  insight  as  to  the  physical  interpretation  of  the  parameter 
6(kT),  a  substitution  of  Eqn  (3-35)  into  Eqn  (3-31)  is  made.  This 
results  in  a  covariance  matrix  update  equation  of  the  form 

P(kT)  =  P{(k-1)T}  -  6(kT)  P{(k-1)T}  T(kT)  TT(kT)  P{(k-1)T}  (3-37) 

From  this  equation,  it  can  be  seen  that  6(kT)  can  be  interpreted  as  a 
gain  term  for  updating  the  covariace  matrix. 

Although  application  of  Eqn  (3-35)  will  result  in  the  desired  value 
of  the  discounting  factor,  because  of  the  restrictions  given  in  Eqn 
(3-33),  Eqn  (3-35)  cannot  always  be  used.  Hagglund  shows  that  by 
incorporating  the  bounds  of  a(kT)  in  Eqn  (3-33)  in  conjunction  with  Eqn 
(3-35)  the  choice  of  a(kT)  becomes 


'(kT)  = 


0 

if 

<  0 

«d(kT) 

if 

0 

<  «j(kT)  ^  -^Yry 

1 

VTkTy 

if 

1 

^klT 

<  «^j(kT)  <  v-l(kT) 

0 

if 

"d 

>  v-l(kT)  + 

(3-38) 


where 


,(kT)  =  TT(kT)  P{(k-1)T}  -.'(kT) 


(3-39) 


The  previous  discussion  is  relevant  to  the  situation  where  the 
parameters  are  constant  or  change  slowly  in  comparison  with  the  time 
constants  of  the  plant.  The  inverse  of  the  covariance  matrix  is  then  a 
good  measure  of  the  information  content.  In  the  case  of  abrupt  plant 
parameter  changes,  P"1(1<T)  is  no  longer  a  good  description  of  the 
information  content  in  the  estimator.  It  will  take  some  time  for  the 
algorithm  to  reflect  the  uncertainty  of  the  old  parameter  estimates  by 
an  increase  of  the  covariance  matrix.  Since  the  covariance  matrix 
partially  determines  the  gain  of  the  of  the  estimation  algorithm,  the 
low  "magnitude"  of  the  P  matrix  will  most  likely  cause  a  slow  parameter 
adaptation  rate.  To  deal  with  this  situation,  Hagglund  developed  a 
"fault"  detection  procedure  to  speed  up  the  adaptation  in  case  of  abrupt 
parameter  changes. 

The  problem  of  how  to  account  for  step-like  parameter  changes  can  be 
broken  down  into  two  parts.  The  first  of  these  is  how  to  detect  the 
parameter  change  (  or  fault  ).  The  second  part  is,  once  the  parameter 
change  has  been  detected,  what  modifications  to  the  estimation  algorithm 
need  to  be  made  to  correctly  account  for  these  changes. 

A  fault  detection  procedure  is  accomplished  by  forming  a  test 
sequence  that  is  sensitive  to  faults.  The  fault  detection  sequence 
should  have  properties  that  are  significantly  different  before  and  after 
the  fault.  Following  the  development  of  a  test  sequence,  the  sequence 
is  evaluated  and  decision  theory  is  applied  to  determine  if  a  fault  has 
occurred. 

Although  the  predominently  used  test  sequence  for  fault  detection  is 


\  -r. 


•  -  «  k  • 


often  the  residual  sequence,  its  use  can  cause  erroneous  results  as 
indicated  in  section  3.3  since  an  increase  in  the  noise  level  can  lead 
to  false  alarms.  To  solve  this  problem,  Hagglund  proposed  to  use  the 
changes  in  the  parameter  estimate  vector  o(kT)  as  the  basis  for  the 
fault  detection  sequence. 

Hagglund  shows  that  when  using  the  differences  between  two  succesive 
^ 

estimates,  AQ(kT)  =  0(kT)  -  G{(k-1)T},  as  the  basis  of  the  fault 

/\ 

detection  sequence,  the  probabilities  of  Ae(kT)  being  positive  versus 
negative  are  approximately  the  same  when  the  estimated  parameters  are 
close  to  the  true  values.  This  assumes  that  the  equation  error  term  in 
Eqns  (3-9),  (3-10),  and  (3-12)  has  white  Gaussian  noise  characteristics. 
Representing  this  mathematically  yields 


P[  A0T(kT)  A0{(k-1)T}  >  0  ]  «  P[  A0T(kT)  A0{(k-1)T}  <  0  ]  (3-40) 

When  a  fault  has  occurred,  the  above  discussion  and  Eqn  (3-40)  no 
longer  hold.  The  estimated  parameters  will  be  driven  toward  the  new 
values  and  therefore  Eqn  (3-40)  will  be  replaced  by 


P[  A0T(kT)  A0{(k-1)T}  >  0  ]  >  P[  A0T(kT)  A0{(k-1)T}  <  0  ]  (3-41) 

Instead  of  observing  the  scalar  product  A0'''(kT)  A0{(k-1)T},  it  is 
more  efficient  to  observe  the  scalar  product  between  A0'^(kT)  and  a  sum 
of  the  latest  estimate  increments.  To  simplify  the  algorithm,  Hagglund 
performs  an  exponential  filtering  of  the  increments  of  the  estimates 
instead  of  an  ordinary  sum.  This  results  in  the  function  w(kT)  being 
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defined  as 


w(kT)  =  >1  w{(k-l)T}  +  A0(kT) 


0  <  yi  <  1  (3-42) 


where  the  design  parameter  -yi  controls  the  number  of  estimate 
increments  that  will  significantly  influence  the  resulting  w(kT).  In 
case  of  a  fault  occurrence,  w(kT)  can  be  viewed  as  being  an  estimate 
of  the  direction  of  the  parameter  change.  This  leads  to  the  test 
sequence  as  developed  by  Hagglund,  s(kT).  The  test  sequence,  s(kT), 
is  defined  to  be 


s(kT)  =  sign[  AQT(kT)  w{(k-l)T}  ]  (3-43) 

The  sign  function  makes  the  test  sequence  insensitive  to  the  noise 
variance.  As  was  the  case  earlier,  when  the  parameter  estimates  are  in 
proximity  to  their  actual  values,  the  function  s(kT)  has  approximately 
a  symmetric  two  point  distribution  with  mass  0.5  each  at  +1  and  -1. 
However,  with  the  occurence  of  a  fault,  the  distribution  is  no  longer 
symmetric  and  the  mass  at  +1  is  larger.  The  idea  behind  the  fault 
detection  technique  is  then  to  inspect  the  latest  values  of  s(kT).  If 
s(kT)  is  +1  an  unlikely  number  of  times  in  a  row,  a  fault  is  then 
declared. 

To  add  the  most  recent  values  of  s(kT),  Hagglund  introduces  the 
stochastic  variable  r(kT)  defined  as 


r(kT)  =  >2  r{(k-l)T}  +  (1-V2)  s(kT) 


0  $  72  <  1  (3-44) 


which  produces  an  exponential  smothing  of  the  test  sequence  s(kT)  in 
order  to  obtain  a  simple  algorithm.  When  the  parameters  are  close  to 
their  true  values,  and  the  equation  error  term  has  "white  noise" 
characteristics,  r(kT)  has  a  mean  value  close  to  zero.  When  a  fault  has 
occurred,  a  positive  mean  is  expected.  A  fault  is  then  declared  when 
the  value  of  r(kT)  exceeds  a  certain  threshold.  The  design  parameter 
here  y2  controls  how  many  s(kT)  values  will  be  included  in  r(kT). 

As  in  the  case  with  most  design  parameters,  there  are  tradeoffs  to 
be  made  when  selecting  the  value  of  >2*  ^  small  value  of  y2  is 
selected,  fast  fault  detection  will  result  at  the  cost  of  less  security 
against  false  alarms.  In  instances  where  the  signal  to  noise  ratio  in 
the  system  is  poor,  speedy  fault  detection  is  not  feasible.  For  cases 
such  as  this,  more  information  is  required  for  proper  decision  making. 
By  assigning  72  ^  larger  value,  more  information  will  be  available  for 
the  fault  detection  algorithm  to  determine  the  presence  of  a  fault. 

As  is  mentioned  earlier,  if  r(kT)  exceeds  a  certain  threshold,  say 
tq,  a  fault  may  be  concluded  with  a  confidence  determined  from  the  value 
of  the  threshold.  The  threshold  rg  can  be  computed  as  a  function  of 
y2  and  the  acceptable  rate  of  false  alarms,  ff,  which  may  be  chosen  to 
suit  the  particular  application  in  question  (19:40).  Hagglund  performed 
this  calculation  and  his  results  are  ilustrated  in  Figure  3-3. 

Having  found  a  way  of  detecting  large  (  and  fast  )  parameter 
changes,  the  next  step  is  to  determine  a  suitable  way  of  increasing  the 
gain  of  the  algorithm  to  speed-up  the  adaptation  rate.  To  accomplish 
this,  Hagglund  chose  to  add  the  quantity  |3(kT)  times  the  identity 
matrix  to  the  covariance  matrix  P(kT).  The  covariance  update  equation 


I 

becomes 


P(kT)  =  P{(I<-1)T}  -  PI(I<  UT)  T(l<TJjfJXkT)  Pi(li:l)T}  ^  j 

[v-l(kT)  -  a(kT)]-l  +  r?(kT) 

where  ?}(kT)  is  given  by  Eqn  (3-39).  The  parameter  /?(kT)  has  the  value 
of  zero  in  all  cases  except  when  a  fault  is  detected.  Then,  a  positive 
value  for  /5(kT)  will  have  the  desired  effect  of  increasing  the 
covariance.  It  remains  to  specify  a  suitable  choice  for  /?(kT). 

One  possibility  of  choosing  /S(kT)  is  to  let  it  be  a  function  of  the 
current  value  of  P(kT)  and  of  the  significance  of  the  fault,  i.e.,  of 
the  value  of  r(kT).  In  deriving  the  equation  for  /?(kT),  Hagglund 
defines  the  following  expression 


/?(kT)  = 


0 


■  v(kT)  . 

T'^CkT)  T(kT) 


[  t..o(kT)  -  i.(kT)  ] 


if  r(kT)  <  rg 

(3-46) 

if  r(kT)  ^  rQ 


where  L.'(kT)  was  shown  to  be  an  eigenvalue  of  the  parameter  update 
equation  that  could  be  assigned  arbitrarily  to  control  the  step  length 
of  the  algorithm  (  a  small  eigenvalue  causes  bigger  steps  towards  the 
new  parameter  values,  while  an  eigenvalue  close  to  one  causes  smaller 
steps  ),  L-o  is  the  current  eigenvalue,  defined  as 

r?(kT) 

^o(kT)  =  1  -  - - -  (3-47) 

v(kT)  +  [  1  -  a(kT)  v(kT)  ]  7?(kT) 


.1 

.1 

.1 
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and  tq  is  the  fault  detection  threshold. 


To  ensure  that  the  covariance  matrix  remains  pcsitive  def inte,  the 
eigenvalue  ^(kT)  needs  to  meet  the  following  condition 

0  ^  ..(kT)  ^  L.o(kT)  (3-43) 

A  simple  way  of  accomplishing  this  is  by  letting  :('<!)  be  a  piecewise 
linear  function  of  the  significance  of  the  fault  as  depicted  in  Figure 
3-4. 


Figure  3-4.  Possible  choice  of  the  eigenvalue  i(kT).  (19:52) 


The  above  choice  of  i'(kT)  causes  Eqn  (3-46)  to  become 


o(kT) 


0 

v(kT)  .qC^T)  (r{(k-l)T}  -  ro) 
YT(kT)  v(kT)  (1  -  rn) 


if  r{(k-l)T}  <  ro 
if  r{(k-l)T}  »  ro 


(3-49) 


I,  '  .  •  »  -*  «•  -  •  .  ,  •  .  •  .  •r,  •  ,  *  ^  •  A  ^  ^  "•'  -,  "■  -«  .  -.  r. 


also  independent  with  respect  to  the  ethers  (  Appendix  C  ),  the 
estimation  problem  can  be  solved  by  adding  the  effects  of  cne  cutpet  at 
a  time  by  going  through  the  algorithm  m  times  per  time  step.  The 
procedure  is  summarized  with  the  following  algorithm 

At  time  kT  (  k  >  0  ,  where  0  is  initiation  time  ) 
calculate  for  i=l,...,m  (see  Eqns  (3-9),  (3-10)  and  (3-12)) 

S'i(kT)  =  S'i.i(kT)  +  Pi(kT)  .i(kT)  ri(kT)  (3-51) 

where 

T-i(kT)  is  the  estimated  prediction  error 
VT(kT)  is  the  estimated  prediction  error  variance 
Pi(kT)  is  the  estimated  parameter  covariance  matrix 
■.■i(kT)  is  the  i^^  column  of  measurements  in  ■r(kT) 

with  initial  conditions 

G'o(O)  initial  presumeo  values  of  the  step-response 
matrix  elements 

v-j(O)  initial  presumed  value  for  prediction  error 
variance  for  i=l , . . .  m 

Pg(0)  estimated  covariance  of  the  par'a”-:-  ' 
at  initiation  time 

■-|(0)  vector  of  past  measurements  prior  to  initiation 
of  identification  for  i=l,...  m 

with  design  parameters 

"a"  desired  variance  of  the  parameter  estimates 

'-'I,  >2»  '"0  design  parameters  for  a  fault  detection  scheme 
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v3  .  design  parameters  of  the  prediction  error 

variance  estimator 


recursive  relationships  (in  proper  order  of  occurrence) 


o(kT) 

=  0'm{(k-l)T} 

(kT) 

=  yi(kT)  -  .iT(kT)  S 

‘i-l(kT) 

l(kT) 

=  Pm{(k-1)T} 

(kT) 

=  LiT(kT)  Pi-i(kT)  . 

i(kT) 

(kT) 

=  vi^kT)  p2i_i(kT) 

ni(kT) 

(kT) 

=  riT(kT)  p3i_i(kT) 

i'i(kT) 

.('■T)  = 


iTkiy 


\i(l<T) 


,.i(kT) 


(kT)  =  rm{(k-l)T} 


(kT)  == 


v-iKk-l)!}  +  (l->3)yi(kT-r)  if  ri_i(kT)  <  r^ 
Vi{(k-1)T}  if  r-i_i(kT)  >  r^ 


,(kT)  =  Vi-l(kT) 


'Hfj(kT) 


''  i  (j  ( k  "i” )  i  ( k  T )  -  1 
if  rtij(kT)«  0 


‘^d 


(kT)  if  0  <  oi^(kT)  «  --y 


'/i 


1 


TTTT  L: 


T:7rirfY  ^  "iri(kT)  <  Vi  ykT) 


(kT) 


I'OiCkT)  =  1  - 


Vi(kT)  +  [  1  -  «i(kT)  Vi(kT)  ]  ,;i(kT) 


(3-63) 


i  \  kT ) 


0  if  r-j-if.kl)  <  rg 

Vi(kT)  .o.(kT)  (ri_i(kT)  -  Rg)  (3-64) 

-  if  r^_i(<T)  j  rg 

vi^(kT)  ri(kT)  (1  -  Rq) 


Pi(kT)  =  Pi.i(kT) 

Pi_l(kT)  .i(kT)  t.iT(kT)  Pi_i(kT) 

- +  .?i(kT)  I  (3-65) 

[vrl(kT)  -  oi(kT)]-l  +  r?i(kT) 

(  equation  (3-51)  evaluated  here  ) 

wo(kT)  =  Wr;,{(k-1)T}  (3-66) 

Wi(kT)  =  VI  Wi.i(kT)  +  [  e'i(kT)  -  Q'i-i(kT)  ]  (3-67) 

Si(<T)  =  sign  [  (  Q'i(kT)  -  G'i.i(kT))T  wi.i(kT)  ]  (3-63) 

ri(kT)  =  >2  ri-i(kT)  +  (1-V2)  Si(kT)  (3-69) 


By  incorporating  Hagglund's  modifications  to  the  recursive 
least-squares  algorithm,  it  is  posible  to  account  for  the  time  variation 
of  plant  parameters  in  a  more  efficient  manner  as  compared  to  the 
standard  RLS  algorithm.  The  major  advantage  of  using  Hagglund's 
algorithm  is  that  it  solves  the  problems  caused  by  nenuniferm  excitation 
of  the  plant.  This  is  important  since  it  allows  for  the  identification 
procedure  to  take  place  during  general  aircraft  maneuvers  by  discounting 
data  depending  on  the  amount  of  information  available.  A  convenient 


feature  is  also  that  the  previcos,  aO  hoc  choices  cf  forgetting  factors 
are  replaced  by  the  perforr.ance  related  choices  cf  the  desired  paraoeter 
variances  (19:115  ). 

3.5  Sc— a-~v 

Chapter  3  has  presented  the  tnecry  tenind  the  reccrsive 
identification  scher.e  that,  together  with  the  control  law  presented  in 
chapter  2,  will  allow  for  the  developn:ent  of  an  adaptive  system  to 
compensate  for  aircraft  parameter  variations.  This  is  done  in  order  to 
demonstrate  the  capability  to  maintain  good  model-following  performance 
while  conducting  in-flight  simulations  under  different  flight 
conditions.  Throughout  these  discussions,  many  variables  have  been 
identified  as  user  defined  design  quantities.  The  next  chapter  will 
provide  the  details  of  the  design  procedures  for  determining  the 
appropriate  design  parameters  fcr  both  the  control  law,  and  the 
recursive  identification  algorithm. 
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IV.  Desi::i  Fr::^'^ 


4. 1  Introd'jction 

The  previous  chapters  in  this  thesis  provioe  sere  insight  into 
the  theory  of  both  the  control  law  and  the  recursive  identification 
algorithm  used  in  the  adaptive  system.  The  intent  of  this  chapter  is 
twofold,  namely,  to  present  a  detailed  description  of  the  procedure  and 
practical  aspects  considered  for  selecting  the  design  parameters  of  the 
adaptive  system,  and  to  provide  details  of  the  system's  representation 
and  simulation  with  MATRIX^. 

Although  deterministic  in  nature,  the  control  law  design  techniques 
developed  by  Prof.  Porter  require  seme  amount  of  trial  and  error  in  the 
development  of  a  particular  design.  A  number  of  variables  are  to  be 
defined  by  the  user  based  on  theoretical  insight  and  experience  gained 
throughout  the  trial  and  error  process.  This  process  must  also  be 
applied  in  adjusting  the  design  parameters  of  Hagglund's  algorithm  to 
optimize  the  adaptation  mechanism's  ability  to  estimate  the  open-loop 
aircraft  dynamics. 

The  adaptive  system  designs  presented  in  this  thesis  are  based  on 
the  use  of  the  linearized  longitudinal  equations  of  motion  of  the 
AFTI/F-16  aircraft  at  a  nominal  flight  condition  of  Mach  0.9  at  10,000 
ft.  of  altitude.  Also,  a  model  of  the  aircraft  at  Mach  0.3,  10,000  ft. 
is  used  to  test  the  system's  ability  to  perform  under  changing 
conditions . 

This  chapter  begins  by  presenting  a  brief  description  of  the 
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"athe”atical  [r^cdels  used  to  rept'esent  tr.e  airsrafp  asd  actuatsr  djras’.ics 
in  section  4.2.  Section  4.3  presents  tne  develcpnents  that  lead  into 
the  final  control  law  design.  General  iratheoatical  considerations,  as 
well  as  peculiarities  of  the  Porter  design  ir.ethcd  are  discussed. 
Finally,  section  4.4  addresses  the  issue  cf  fine  tuning  the  estimator 
design  variables,  and  also  highlights  the  various  practical  signal 
processing  aspects  required  for  its  implementation. 

4.2  Aircraft  Model 

The  airplane  dynamics  used  in  this  study  are  represented  by  a 
set  of  first  order  matrix  diferential  equations  in  the  state  space  form 
of  equations  (2-1)  and  (2-9).  These  state  space  equations  are  obtained 
from  the  forces  and  moments  acting  upon  the  aircraft,  and  are  expressed 
using  the  body  axis  system  centered  at  the  aircraft's  center  of  gravity 
(Figure  4-1).  A  detailed  derivation  of  the  longitudinal  state 
perturbation  equations  used  in  this  thesis  is  presented  in  Reference  7. 

Equation  (4-1)  shows  the  longitudinal  state  space  mcdel  c*  r  - 
AFTI/F-16  in  terms  of  the  primed  dimensional  iced  staoility 


e  IS  the  pitch  angle 

u  is  the  forward  velocity 

a  is  the  angle  of  attack 

q  is  the  pitch  rate 

6q  is  the  elevator  deflection 
6f  is  the  flaperon  deflection 


The  equations  used  to  calculate  the  coefficients  in  Eqn  (4-1)  are  also 
found  in  Reference  7.  In  this  case,  a  Flight  Dynamics  Laboratory 
aerodynamic  data  package  for  the  AFTI/F-16  was  used  to  obtain  the 
coeficients  in  Eqn  (4-1)  by  trimming  the  aircraft  at  the  previously 
named  flight  conditions.  The  relevant  data  for  these  models  is 
presented  in  appendix  A. 

The  next  step  is  simply  to  define  the  output  vector  of  the  quantites 
of  interest.  In  this  case,  the  desired  outputs  are  flight  path  angle 
and  pitch  rate.  The  flight  path  angle  is  defined  as 


(4-2) 


This  leads  to  the  following  output  equation 


These  equations,  representing  the  aircraft  dynamics,  are  implemented 
by  means  of  the  computer  aided  design  package  MATRIX^  with  its 
simulation  facility  SYSTEM-BUILD  (23).  This  is  ilustrated  in  Figure  4-2 
thru  4-4.  Figure  4-2  shows  the  elements  of  the  block  "A/C"  wich  depicts 
the  overall  aircraft  model  with  a  representation  directly  derived  from 
the  state  space  equations.  The  different  set  of  dynamics  associated 
with  the  two  flight  conditons  of  interest  are  implemented  by  means  of 
gain  scheduler  functions  (Figures  4-3  -  4-4)  which  change  the  elements 
of  the  A  and  B  matrices  in  the  model  of  Figure  4-2.  The  elements  in 
the  C  matrix  remain  constant  throughout  the  simulation. 

Figure  4-2  also  shows  the  addition  of  actuator  dynamics  into  the 
aircraft  simulation.  These  include  nonlinearities  such  as  surface 
position  and  rate  limits.  Consideration  of  their  effects  is  very 
important,  especially  in  the  case  were  the  control  system  is  to  provide 
artificial  stability  to  the  aircraft  as  it  is  the  situation  here. 

The  actuator  model  used  in  this  simulation  is  derived  from  the  VISTA 
design  requirements  that  call  for  faster  actuator  dynamics  than  those 
found  in  the  regular  F-16  aircraft.  For  purposes  of  these  simulations  a 
simple  first  order  model  will  be  used.  The  actuator  model  for  both  the 
elevator  and  flaperon  is  of  the  form 

^  _  44  _ 

‘■•cmd  44ir 

or  equivalently 

,s  =  -44  0  +  44 


(4-4a) 

(4-4b) 


Equation  (4-4)  is  implemented  for  both  actuators  with  the  models  shown 


(  C  o  n  T  I  n  u  a  iJ  *i  ) 


Figure  4-2.  Aircraft  Dynamics  Representation 
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Figure  4-4.  Gain  Scheduler  Function  for  the  Aircraft's  B  Matrix 


derived  from  Eqn  (4-4b),  with  the  addition  of  the  rate  and  position 
limiting  functions. 

The  control  surface  rate  is  limited  simply  by  placing  a  limiter 
function  after  the  rate  of  Eqn  (4-4b)  is  calculated.  This  is 
implemented  after  the  summing  junction  in  Figure  4-5.  The 
implementation  of  the  control  surface  position  limiting  function  is 
slightly  more  involved.  It  would  be  erroneous  to  simply  place  a  limiter 
after  the  position  of  the  surface  is  calculated  from  the  integration  of 
the  surface  rate,  since  this  would  imply  that  the  calculated  surface 
rate  of  motion  would  not  be  zero  when  the  surface  is  stuck  at  a  limit. 
Instead,  a  limited  integrator  function  is  implemented  to  calculate  the 
control  surface  position  from  the  calculated  surface  rate.  This  is 
ilustrated  in  Figure  4-6.  The  rate  canceler  in  the  limiter  integrator 
(block  labeled  R/CA'JC)  monitors  the  calculated  surface  position  and 
compares  it  to  the  surface's  position  limits.  If  the  surface  deflection 
is  at  a  limit,  the  rate  canceler  function  negates  the  commanded  rate  at 
the  summing  junction  location.  This  causes  the  desired  effects  of 
zeroing  the  surface  rate  and  stoping  the  surface  position  calculation  at 
its  current  point  until  the  commanded  surface  rate  changes  direction. 
While  the  surface  deflections  are  within  the  normal  range,  the  rale 
canceler  simply  outputs  zero,  thus  letting  the  incoming  surface  rate  be 
integrated  to  obtain  the  surface  position.  The  control  surface  position 
and  rate  limits  used  in  this  study  are  given  in  Table  4-1. 

The  position  limits  are  relative  from  the  trim  position  of  the 
surfaces.  Since  the  simulation  is  based  on  perturbation  equation  of 
motion  at  different  flight  conditions,  these  limits  are  scheduled 
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Table  4-1 


Control  Surface  Position  and  Rate  Limits 


Surface 

Fit.  Condition 

Position  Limit 
(deg) 

Rate  Limit 
(deg/sec) 

Elevator 

Nominal 

+  27.37/  -  22.63 

90 

II 

Off  Nominal 

+  27.06/  -  22.94 

II 

Flaperon 

Nominal 

+  22.0  /  -  21.0 

78 

II 

Off  nominal 

+  7.54/  -  35.46 

II 

according  to  which  model  is  being  used  at  any  given  time.  The  absolute 
limits  are  *  25  deg.  for  the  elevator  and  +  20  deg.  /  -  23  deg.  for  the 
flaperon.  The  direction  of  positive  surface  deflection  Is  indicated  in 
Figure  4-7. 

Finally,  consideration  to  imperfect  measurements  is  given  in  the 
design  by  providing  the  capability  to  inject  measurement  noise  into  the 
simulation.  This  is  ilustrated  in  Figure  4-8,  which  shows  the  details 
of  the  sensor  block.  A  simple  model  for  zero-mean,  white,  Gaussian 
noise  is  used  to  corrupt  the  individual  measurements  independently.  The 
specific  noise  levels  used  are  provided  in  the  next  chapter.  Figure  4-8 
also  shows  the  capability  of  including  sensor  dynamics  in  the 
simulation.  The  sensor  dynamics  implementation  in  this  case  is  limited 
to  simple  low-pass  anti-aliasing  filters  with  corner  frequency  of  40  Hz. 
The  reason  for  this  is  that,  at  the  time  of  this  study,  no  specific 
information  was  available  on  the  sensors  to  be  used  in  the  design  of  the 
variab''e  stability  system  for  VISTA. 
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Figure  4-8.  Sensor  Measurement  Noise  Representation 
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4.3  Control  Law  Design 

4.3.1  Mathematical  Background.  Before  a  particular  control  law 
design  is  developed  it  is  necessary  to  make  sure  that  the  system  for 
which  it  is  intended  satisfies  a  number  of  requirements.  Of  fundamental 
importance  is  for  the  plant  to  be  completely  controllable  and 
observable.  Controllability  implies  that  the  inputs  can  affect  all  the 
modes  (poles)  of  the  system,  thus  making  it  possible  to  alter  these 
modes  in  order  to  effectively  modify  the  system's  characteristics.  For 
a  linear  time-invariant  system  described  by  Eqn  (2-1)  controllability 
can  be  checked  by  evaluating  the  rank  of  the  controllability  matrix 
(14,40),  denoted  M^-,  given  by 


Mr  =  [  B  I  AB  I  a2b  I  ...  I  A^-Ib  ] 


(4-5) 


where  n  is  the  number  of  states.  If  the  rank  of  Mq  is  equal  to  n, 
then  the  system  is  completely  controllable.  If  M(-  is  rank  defective  the 
system  is  deemed  uncontrollable  signifying  that  some  of  the  modes  cannot 
be  affected  by  the  input  to  the  plant.  In  this  case  the  rank  defect  of 
M(-  (i.e.,  n  minus  the  actual  rank  of  Mq)  tells  us  how  many  modes  are 
uncontrolable. 

Observability  implies  that  the  outputs  of  the  plant  are  affected  by 
e;e’'y  mode,  and  that  every  state  affect  the  outputs  of  the  plant  (40). 
This  property  is  particularly  important  since  control  is  to  be  provided 
by  cutp'jt  feedback.  Although  not  all  the  states  are  to  be  measured, 
their  influence  on  the  outputs  is  necessary  if  proper  control  of  the 
modes  is  to  be  achieved.  In  addition,  observability  becomes  necessity 
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in  this  particular  case  where  the  dynamics  of  the  plant  are  to  be 
identified  from  input-output  data  for  the  purpose  of  applying  adaptive 
control.  Lack  of  observability  would  imply  that  the  plant  could  not  be 
properly  characterized  solely  from  input-output  data.  This  would  lead 
to  an  erroneous  plant  representation  and,  in  turn,  an  incorrect  control 
strategy.  For  the  linear  time-invariant  system  represented  by  Eqns 
(2-1)  and  (2-9)  the  property  of  observability  can  be  checked  by 
evaluating  the  rank  of  the  observability  matrix  Mq  (14,40),  defined  as 

Mq  =  [  CT  I  aTcT  I  (aT)2cT  i  ,  (AT)n-lcT  ]  (4-6) 

Analogously  to  the  controllability  check,  if  the  rank  of  Mq  is  equal 
to  n  then  the  plant  is  deemed  fully  observable.  If  this  is  not  the 
case,  then  the  rank  defect  of  Mq  indicates  the  number  of  unobservable 
modes . 

Controllability  and  observability  are  checked  for  the  two  sets  of 
dynamics  used  in  this  thesis  (Appendix  A).  This  is  done  with  the  aid  of 
the  matrix  manipulation  and  control  design  functions  of  the  MATRIX^ 
package.  Both  models  satisfy  the  full  controllability  and  observability 
conditions . 

Another  important  consideration  in  the  Porter  method  is  that  of  the 
location  of  transmission  zeros.  The  transmission  zeros  of  the  plant, 
defined  as  zeros  of  the  equivalent  transfer  function  representation  that 
block  particular  modes  frcm  the  input,  are  regions  to  which  some  of  the 
slow  roots  of  the  closed  loop  system  migrate  as  the  gain  approaches 


infinity.  Output  feedback  does  not  alter  the  location  of  transmission 


vs 


zeros,  therefore  it  is  desired  that  these  zeros  be  stable  in  order  to 
insure  system  stability  at  high  values  of  gain.  This  however  does  not 
guarantee  stability  for  lower  values  of  gain  for  which  the  system  roots 
may  pass  through  unstable  regions  before  approaching  the  transmission 
zero  locations.  It  is  therefore  imperative  that  the  system's  roots  be 
checked  throughout  the  design  stages  to  verify  the  stable  location  of 
the  closed-loop  roots.  This  task  is  easily  accomplished  with  the 
analysis  tools  provided  by  MATRIXy. 

Of  particular  interest  is  the  existence  of  a  transmission  zero  at 
the  origin.  This  condition  results  from  the  inclusion  of  the  pitch  rate 
in  the  output  vector.  The  presence  of  a  transmission  zero  at  the  origin 
suggest  that  the  system  is  unstable  and  uncontrollable  (35)  because  of 
the  location  of  a  closed-loop  system  root  at  the  transmission  zero 
location.  An  analysis  of  this  situation  is  required  at  this  point.  As 
pointed  out  by  Barfield  (7:83-85)  the  transmission  zero  at  the  origin 
results  from  the  assumptions  made  in  forming  the  equations  of  motion  for 
the  aircraft  in  which  the  following  relationship  between  the  pitch  rate 
and  pitch  angle  is  established 


(4-7) 


If  q  is  to  be  commanded  with  step  functions,  would  ramp  to  infinity 
due  to  the  integration  caused  by  the  pole  located  at  the  transmission 
zero  location.  Thus  seme  of  the  system's  responses  are  unbounded  for 
bounded  inputs.  For  any  practical  maneuver,  q  cen—ands  resemble  pulses, 
more  than  steps,  or  perhaps  some  other  shape  applied  for  a  finite  amount 
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of  tine  that  produces  a  bounded  pitch  angle  response.  This  is  precisely 
what  is  expected  given  the  relationship  of  Eqn  (125). 

A  by-product  resulting  from  the  eA-istance  of  the  transmission  zero 
at  the  origin  is  that  the  introduction  of  integral  action  in  the  control 
law  produces  what  is  termed  a  "functionally  uncontrollable"  system  (35). 
The  vector  integrator  that  is  part  of  the  cont^'ol  law  introduces 
additional  dynamics  in  the  control  loop.  As  part  of  the  design 
procedure,  it  is  necessary  to  check  for  controllability  of  the  augmented 
system  composed  of  the  control  law  and  plant  dynamics  (which  may  include 
the  actuators  and  sensors).  A  simple  check  for  this  is  performed  by 
satisfying  the  following  rank  test  (38) 


:  =  n  +  m 


(4-8) 


Executing  this  test  for  the  aircraft  moaels  of  this  thesis  results  in  a 
rank  deficiency  of  one  for  both  models  in  question.  Again,  additional 
discussion  is  required  at  this  tim.e.  A  rank  deficiency  of  one  implies 
that  only  one  closed-loop  pole  is  uncontrollable.  An  analysis  of  the 
behavior  of  the  closed-loop  poles  of  the  system  (Appendix  D)  reveals 
that  the  uncontrol  Table  mode  corresponds  to  one  of  the  poles  at  the 
origin  introduced  by  the  vector  integrator  of  the  control  law.  This 
uncontrollable  closed-loop  pole  at  the  origin  accounts  for  the 
integration  that  takes  place  due  to  Eqn  (4-7).  The  fact  that  the 
location  of  this  root  cannot  be  changed  merely  reflects  the  case  that 
the  definition  given  in  Eqn  (4-7)  must  prevail  in  the  resulting  pitch 


angle  response.  The  other  roots  of  the  system  are  controllable  and 
their  asymptotic  behavior  follo'^s  the  theory  presented  in  chapter  2. 
The  condition  of  functional  uncontrollability  is  not  considered  an 
adverse  indication  in  this  design  but  rather  a  reflection  of  the 
assumptions  made  in  defining  the  dynamics  of  the  problem.  The  system  is 
then  deemed  conditonally  stable  for  the  particular  design  in  question, 
and  thus  considered  satisfactory  for  purposes  of  this  study. 

4.3.2  Design  Variables.  The  key  elements  of  the  proportional 
plus  integral  control  law  of  Eqn  (2-46)  are  the  proportional  gain  matrix 
K,  and  the  integral  gain  matrix  K  .  Once  the  step-response  matrix  of 
the  aircraft  model  is  known  for  a  particular  flight  condition  and 
sampling  time,  the  gain  matrices  are  calculated  according  to  Eqns  (2-47) 
and  (2-48)  using  the  following  design  variables: 

:  =  the  diagonal  weighting  matrix,  diag{  r  ,  r  } 

.  =  the  ratio  of  integral  to  proportional  control 

The  values  of  c  ,  c  ,  and  .  are  chosen  so  that  se.eral  cbjectives  are 
satisfied;  the  first  of  which  is  to  stabilise  the  aircraft.  After  this 
is  accomplished  the  controller  is  tailored  to  produce  the  desired 
tracking  characteristics.  The  selection  c*  tn.ece  c-.-sign  par'a-.ete’'S  in 
this  thesis  is  based  on  the  use  cf  an  aircj'aft  mcd-'  tri'''~-c  at  Mach  0.9 
and  10,0G0  ft.  altitude,  with  a  sa-^pling  time  cf  0.01  sec.  for  the 
control  law.  The  value  chosen  for  the  sa—pli":  t'"'e  is  conside''ed 
sufficiently  snort  as  to  be  p;ssipie  to  oos-'-.e  a’'  re>;-.a-':  moCes  cf 
the  open-loop  plant,  a-od  is  w'thin  t:'e  capap-^l'ty  of  ti'e 
state-cf-the-art  in  flimot  cc'.ti'ol  co"'p-t^’'S. 


The  Porter  design  technique  is  based  primarily  on  the  assumption  of 
linear  plant  models,  and  the  availability  of  sufficient  control 
authority.  The  limitations  imposed  by  control  surface  deflection  and 


rate  limits  make  the  system  nonlinear.  This  situation  imposes 
additional  constraints  in  the  design  task.  The  problem  then  becomes  the 
formulation  of  a  control  law  that  provides  for  sufficiently  fast 
responses,  with  an  acceptable  level  of  tracking  error,  and  while  staying 
within  control  surface  deflection  and  rate  limits.  A  trial  and  error 
procedure  is  used  to  achieve  this  purpose. 

Initial  guidance  in  the  trial  and  error  process  can  be  gained  by 
considering  the  asymptotic  characteristics  of  the  closed-loop  system. 
Recalling  from  the  discussion  in  chapter  2,  as  the  sampling  frequncy  is 
increased  the  closed-loop  tranfer  function  matrix  (Eqn  (2-42)) 
approaches  the  diagonal  form 


0 
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r(z)  =  diag  ■,  :  (*^-9) 

which  contains  only  the  fast  roots  of  the  system.  This  relationship 
provides  the  insight  as  to  how  the  system  exhibits  increasingly  tight 
and  decoupled  control  witn  inc^'eas ing  gain.  It  is  cbvious  then  that  by 
proper  selection  of  the  elements  of  the  t  weighting  matrix  the 
transient  responses  can  be  independently  adjusted,  and  thus  the  speed  at 
which  the  cutouts  follow  f'e’r  corresponding  cc"'-and  signal  from  the 
model  aircra-'t.  Tne  leoel  cf  e^'ror  in  tne  responses  can  then  be 
adjusted  b>  the  se'ccticn  of  tne  ratio  cf  c^'cp:''ticnal  to  integral 
control  (  .  ).  This  rati’.s  also  cc-'t-ols  s:"’e  c"  f  e  slew  roots  cf  the 


system.  Although  not  present  in  the  asymptotic  transfer  function,  at 
finite  values  of  gain  these  slow  roots  become  observable  and  thereby  are 
present  in  the  outputs. 

The  control  law  representation  in  MATRIX^  is  given  in  Figure  4-9. 
Its  implementation  follows  directly  from  Eqn  (2-46)  with  the  exception 
of  the  block  labeled  "INTLIM"  that  stands  for  integration  limiter.  The 
intent  of  this  function  is  to  stop  the  integration  of  the  error  signals 
if  any  control  surface  is  commanded  to  a  position  limit.  The  integration 
of  the  error  signals  continue  once  the  affected  surface(s)  leave  the 
maximun  deflection  position.  This  action  can  improve  the  closed-loop 
system  response  signicantly  (see  Appendix  D)  in  a  manner  analogous  to 
that  discussed  in  the  implementation  of  the  actuator  models. 

All  of  the  elements  of  the  MATRIX^  implemetation  discussed  so  far 
are  now  collected  in  the  closed-loop  simulation  diagram  of  Figure  4-10. 
The  block  labeled  AOAPTZ  implements  the  recursive  identification 
algorithm  and  is  discussed  in  section  4.4.  The  simulation  setup  of 
Figure  4-10  permits  a  trial  and  error  evaluation  of  different 
combinations  of  weighting  matrix  elements  and  proportional  to  integral 
control  ratios.  At  each  design  iteration,  and  prior  to  any  simulation, 
the  poles  of  the  closed-loop  system  are  checked  as  a  first  means  of 
assesing  the  the  stability  of  the  design.  If  the  particular  selection 
of  design  parameters  produces  a  stable  location  for  the  closed-loop 
poles,  then  a  time  response  analysis  is  conducted  via  simulation.  This 
is  necessary  since  tne  close-locp  pole  location  analysis  can  only  be 
dene  for  the  linear  representation  of  the  sgstem.  The  nonlinear  effects 
can  be  better  acounted  for  in  simulations. 


4.3.3  Command  Maneuver.  In  the  case  of  a  model-following 
control  law,  the  command  inputs  are  typically  a  set  of  time  histories 
coming  from  a  model  of  the  simulated  vehicle  that  is  being  tested  by  the 
evaluation  pilot.  These  model  responses  have  to  be  modified  by  a  number 
of  transformations  that  account  for  differences  in  geometry,  attitude, 
and  speed  between  the  model  and  host  aircraft  (Figure  4-11).  To 
simplify  the  design  procedure,  this  thesis  addresses  the  case  in  which 
the  host  aircraft  is  following  a  model  of  itself  at  an  identical  flight 
condition  throughout  a  particular  maneuver.  This  has  the  effect  of 
making  the  aforementioned  transformation  equal  to  identity.  This  action 
is  practically  motivated  since,  in  the  case  of  the  actual  flight  test 
program  of  the  VISTA  it  will  provide  a  highly  controlled  experiment  to 
evaluate  and  validate  the  model  of  the  host  vehicle.  The  model  of  the 
host  vehicle  will  be  required  for  ground  simulation  efforts  and  variable 
stability  system  control  law  development. 

The  flight  path  angle  and  pitch  rate  time  histories  used  in  this 
thesis  were  obtained  from  a  real-time  simulation  of  the  AFTI/F-16, 
conducted  at  the  Flight  Dynamics  Laboratory's  Control  Synthesis  Branch. 
These  time  histories  are  shown  in  Figures  4-12  and  4-13. 
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Figures  4-14  thru  4-69  show  a  representative  sample  of  the  design 
simulations  conducted. 

To  assess  the  tracking  performance  of  the  system,  a  model  following 
performance  index  is  calculated  for  each  of  the  commanded  quantities. 
It  is  desired  that  the  average  response  error  absolute  value  of  each 
output  quantity  should  not  exceed  the  average  absolute  value  of  ten 
percent  of  its  corresponding  reference  signal.  This  is  expressed  as 
fol lows 


I 

yni(T)  -  yh(i)  I  I  .l(yni(i))  |  dt  (4-10) 


where  the  subscript  m  indicates  a  model  output,  the  subscript  h 
indicates  a  host  airplane  output,  tf  the  ending  time  of  the  maneuver, 
and  the  index  i  pertaints  to  the  i^^  element  of  the  output  vector.  In 
addition  to  a  model-following  tracking  criteria,  the  control  surface 
deflections  are  checked  for  possible  reaching  of  position  and/or  rate 
limits. 

Based  on  the  results  of  these  simulations  a  final  controller  design 
is  chosen  that  satisfies  the  desired  performance  specifications  without 
putting  serious  demands  upon  the  actuator  deflections  and  rates  for  the 
specific  maneuver  used.  The  selected  controller  design  parameters  are 


given  in  Table  4-2. 


Table  4-2 


Final  Control  Law  Design  Parameters 


I  0.3 


=  0.8 


Figure 


-14.  Flight  path  angle  command  and  response  (deg) 
I  =  diag{  0.1,  0.7  } ,  p  =  0.8 


Fiaure  4-15.  Flight  path  angle  tracking  performance  index  ( 
r  =  diag{  0.1,  0.7  },  .■  =  0.8 
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Figure  4-28.  Flaperon  deflection  (deg), 
:  =  diag{  0.3,  0.7  },  ,■  =  0.8 
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Figure  4-29.  Flaperon  deflection  rate  (deg/sec), 
:  =  diag{  0.3,  0.7  } ,  .  =0.8 
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Figure  4-42.  Elevator  deflection  (deg). 
L  =  diag{  0.3,  0.5  },  p  =  0.8 


4-43.  Elevator  deflection  rate  (ceg/sec). 
L  =  dia-f  0.3,  0.5  },  ,  =  0.3 
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Figure  4-44.  Flaperon  deflection  (deg). 
I  =  diag{  0.3,  0.5  },  r  =  0.8 


Figure  4-45.  Flapercn  deflection  rate  (deg/sec), 
:  =  diag(  0.3,  0.5  },  .  =  O.S 


Figure  4-48.  Pitch  rate  coirrrand  and  response  (deg/sec). 
L  =  diag{  0.3,  0.8  } ,  <  =  0.8 
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Figure  4-49.  Pitch  rate  tracking  perforn^ance  index  (deg/sec) 
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4.3.5  Parameter  Variation  and  Sensor  Noise.  To  assess  the 

capability  of  the  resulting  fixed  gain  control  law  design  in  the 

presence  of  plant  uncertainty,  the  controller  resulting  frcm  the  design 
parameters  of  Table  4-2  is  subjected  to  step  changes  in  the  plant 
dynamics  by  using  a  model  of  the  AFTI/F-16  trimmed  at  a  flight  condition 
of  rrach  0.3  and  10,000  ft.  MSL.  This  is  to  represent  either  uncertainty 
on  the  knowledge  of  the  stability  derivatives  of  the  aircraft,  or  sim.ply 
a  change  of  flight  condition  of  the  host  aircraft  during  a  flight.  In 

addition,  several  simulations  are  conducted  with  the  inclusion  of 

various  levels  of  sensor  measurement  noise.  This  is  done  either 
independently,  or  jointly  with  plant  step  changes  in  order  to  determine 
their  effect  on  simulation  fidelity.  Results  of  these  tests  are 
d:cu"ented  in  chapter  5. 

. 4  F Ice-'t "'f ie*'  Desicn 

4.4.1  Cesic"  Va’~iables.  Although  seemingly  complex,  the 
!'t.cursi-e  identif icaticn  algorithm  proposed  by  Hagglund  has  relatively 
fcw  design  variacles.  Their  specification  is  mostly  based  on  experience 
gain.ed  cy  trial  and  errcr  and  a  few  analytical  relationships  as  it  is 
si'iOAn  ncAt.  Recalling  from  chapter  3,  these  design  param^eters  are 

desired  variance  cf  the  parameter  estimates 

design  para~eters  for  a  fault  detection  sche'”e 

design  paramecers  of  the  pt'edicticn  error 
variance  estimatcr 
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The  next  paragraphs  will  ad 
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nt  aspects  involved 
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in  their  final  selection. 


4.4.2  Parai~eter  Estirrate  Variance.  The  first  of  the  design 
parameters  to  be  considered  is  the  desired  variance  of  the  parameter 
estim.ates  (  "a"  in  Eq  (3-53)  ).  Simply  stated,  this  design  parameter 
correspond  to  a  specification  of  how  precise  we  wish  the  parameter 
estimates  co  be.  The  selection  of  "a"  is  conceptually  a  simpler 
choice  cc  , pared  to  the  selection  of  the  forgetting  factor  in  the 
standard  RLS  algorihtm,  since  it  is  not  based  on  any  assunptions  of  the 
rate  of  change  of  plant  parameters. 

Typical  considerations  for  the  selection  of  "a"  are  the  relative 
magnitudes  of  the  parameters  to  be  estimated  and  the  noise  level  in  the 
input-output  data  used  in  the  algorithm.  It  may  be  argued  that  these 
considerations  are  themselves  subject  of  speculation  by  the  designer  in 
the  case  of  unknown  plants,  howeve'*,  this  issue  is  not  ccnsidered  a 
critical  one  due  to  a  peculiar  property  of  the  algcrihcm  used.  Given  a 
particular  level  of  plant  eAcitaticn  and  signal  to  noise  ratio,  the 
algorithm  does  its  best  to  converge  to  tl'e  specified  parameter  variance 
by  adjusting  its  time  horiccn  for  convergence.  Stability  analysis 
conducted  by  Hagglund  indicate  however  that  if  the  signal  to  noise  ratio 
happens  to  be  lar'ger  than  the  one  assumed,  the  F  matrix  will  net 
convenge  to  a-I,  but  it  will  be  actually  smaller  (19:61).  This  is 
hardly  any  problem  since  the  algorithm  will  then  pT'ovice  a  higher 
accuracy  for  the  estimates  than  the  one  specified.  A  few  design 
iterations  can  easily  reveal  a  suitable  choice  fer  "a\  The  particular 
value  selected  for  this  investigation  is  5.0  E-5. 


4.4.3  Fault  Detector  Design.  The  next  design  step  is  the 
selection  of  suitable  parameters  for  the  fault  detection  algorithm.  Its 
design  begins  with  the  selection  of  the  integrator  gain  in  Eqns 
(3-42)  and  (3-67).  The  value  of  >1  controls  the  number  of  most  recent 
parameter  updates  that  affect  the  estimation  of  the  average  direction  in 
which  the  parameter  vector  is  updated.  It  is  desired  for  to  be 
relatively  high  (  close  to  1  )  so  that  a  sufficient  number  of  samples 
are  included  in  the  estimate  of  the  update  direction,  but  not  so  high 
that  it  would  delay  the  detection  of  direction  change  due  to  plant 
parameter  changes.  After  conducting  various  simulations,  the  value  of 
•yi  =  0.85  was  found  to  give  satisfactory  results  for  the  particular 
case  considered  here. 

The  selection  of  >2  (3-44)  and  (3-69)  is  based  on  the  same 
rationale  used  in  selecting  yj.  Hagglund  indicates  (19:39)  that  a 
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the  most  recent  values  of  s(kT)  in  determining  the  stochastic  variable 
r(kT).  This  is  considered  satisfactory  for  this  study  since  it  only 
introduces  a  detection  delay  between  20  and  30  samples. 

The  remaining  design  variable,  the  threshold  rg,  is  assigned  a 
value  that  depends  on  the  value  chosen  for  )/2»  what  is  considered 
an  acceptable  rate  of  false  alarms  for  the  particular  design  in 
question.  If  a  small  value  is  selected  for  the  threshold  it  is  possible 
to  detect  faults  quickly,  but  the  false  alarm  rate  will  be  high.  The 
chart  produced  by  Hagglund,  ilustrated  in  Figure  3-3,  shows  the 
compromise  on  false  alarm  rate  versus  the  value  of  the  threshold  rQ, 
for  various  values  of  >2.  For  this  study  it  is  considered  that  an 
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expected  false  alarm  rate  of  one  in  one  thousand  samples  is 
satisfactory,  given  the  length  of  the  design  simulations.  Using  Figure 
3-3,  and  the  values  chosen  for  >2  and  ff,  a  value  of  rg  =  0.5  is 
selected. 

4.4.4  Noise  variance  estimator  design.  The  design  of  the 
noise  variance  estimator  is  composed  of  the  selection  of  the  residual 
weighting  coeficient  >3,  the  time  delay  compensation  parameter  7,  and 
the  threshold  r^  in  Eqns  (3-50)  and  (3-60).  Initially,  Eqn  (3-50)  is 
evaluated  independently  in  several  design  simulations  with  a  Gaussian 
white  noise  source  of  known  variance.  These  simulations  are  used  to 
determine  a  suitable  value  for  >3  that  produces  results  as  close  as 
posible  to  the  known  variance  of  the  noise  generator.  This  is  done  with 
7  set  to  zero.  A  value  of  >3  =  0.95  is  found  to  produce  satisfactory 
results. 

The  time  delay  compensation  parameter  r  is  easily  adjusted, 
recalling  that  the  selection  >2  previously  mentioned  caused  a 
detection  delay  of  approximately  20  to  30  samples.  The  parameter  7  is 
therefore  assigned  a  value  of  20.  Finally,  it  remains  to  assign  a  value 
for  the  threshold  rj^.  This  is  done  on  a  trial  an  error  basis  by 
repeated  simulation  with  the  entire  algorithm,  a"d  noticing  how  the 
f luct'..ations  in  the  parameter  estim.ates  affect  the  signal  r(kT).  The 
threshold  rj^  is  usually  chosen  smaller'  than  rg  since  it  is  desired  to 
exclude  the  effects  of  parameter  estimate  tt'ansienps  on  the  residuals 
used  to  estimate  the  noise  variance.  These  transients  are  typical  in 
this  type  of  estimation  algorithms,  but  are  not  necessarily  the  result 
of  plant  parameter  changes.  The  threshold  ri  =  0.2  is  used  in  this 


Trns  concludes  the  dicussion  on  the  selection  of  the  primary  design 


variables  of  the  identification  algorithm.  Several  practical 
considerations  however,  make  necessary  the  use  of  additional  signal 
processing  techniques  to  complete  the  implementation  of  the 
identification  algorithm.  Section  4.4.5  discusses  this  topic. 

4.4.5  Practical  Signal  Processing  Considerations.  In  order  to 
implement  (  and  in  some  instances  improve  the  performance  of  )  the 
identification  algorithm,  it  is  necesary  to  use  a  series  of  additional 
filters  to  ensure  that  the  signals  used  in  the  identification  and  the 
control  law  design  process  are  well  conditioned  for  such  purposes.  For 
example,  the  estimation  algorithm  is  designed  to  use  perturbations  of 
the  plant's  input-output  signals,  however,  the  sensors  typically  provide 
data  on  "absolute"  measurements.  In  other  words,  the  measurements 
U(kT)  and  Y(kT)  are  composed  of  nominal  values  U.>(kT)  and  y,.(kT), 
and  perturbation  signals  u(kT)  and  y(kT).  This  can  be  expressed  as 


U(kT)  =  U.,(kT)  +  t;(kT)  (4-11) 

Y(kT)  =  Y.(kT)  +  y(kT)  (4-12) 


In  the  case  cf  the 
as3./-'ea  as  tne  coma 
obtain  the  pe’'turba 
U , (kT)  would  involve 
what  w6  ar'e  tryini 


suto'-t  signal,  the  ncs'inal  value  Y.,(kT)  may  be 
inded  quantity.  Then  Eqn  (4-12)  may  be  used  to 
icns.  Unfotur.ately ,  tne  prccedure  for  obtaining 
inverting  the  model  of  the  plant  which  is  of  course 
to  determine.  It  is  hianl.  desit'able  to  use  a 


technique  that  does  not  requires  knowledge  of  the  nominal  measurements. 
Assuming  that  any  changes  in  the  nominal  values  (caused  by  the  maneuvers 
of  the  aircraft)  are  slow  compared  to  the  sampling  period,  so  that  they 
remain  relatively  constant  between  any  two  consecutive  sampling 
instants,  the  perturbations  may  be  approximated  using  the  differences 


U(kT)  -  U{(k-1)T]  «  u{kJ)  -  u{{k-l)J}  «  ^uikJ) 


(4-13) 


Y(kT)  -  Y{(k-1)T}  -  y(kT)  -  y{(k-l)T}  =  Ay(kT) 


(4-14) 


These  differences  essentially  produce  a  high-pass  filtering  effect  on 
the  measured  quantities,  thus  removing  their  nominal  components. 
Instead  of  u{kT)  and  y(kT),  the  signals  iu{kl)  and  Ay(kT)  are 
used  for  the  parameter  estimation  (25). 

Another  consideration  in  the  ccnciticning  of  input-output  data, 
especially  in  the  case  of  a  tracking  control  law,  is  the  effect  that 
changes  in  the  command  signal  have  on  the  parameter  estimates.  Because 
of  the  high-pass  filtering  effect  of  Bqns  (4-13)  and  (4-14),  abrupt 
changes  in  the  command  signals  to  the  closed-loop  system  produce  spikes 
in  the  input-output  differences  Au(kT)  and  ,\y(kT).  These  spikes  are 
reflected  in  the  residuals  that  are  used  to  update  the  par'a"'eter  and 
noise  variance  estimates  by  virtue  of  Eqns  (3-51),  (3-53),  and  (3-60); 
thus  producing  spikes  in  the  para-eter  estimates  the-^sel ves .  This 
effect  is  highly  undesirable,  however,  a  reduction  of  t’:e  para~eter 
variations  can  be  achieved  by  filtering  both  the  inpot  and  output 
differences  with  identical  lew  pass  filter  al  gcr'i  t:'"'S .  Vat'icus  forms  of 


lew  pass  aldcrithms  may  be  used  for  this 


"e  c-e  selected  for' 


this  application  is  of  the  form 


fcvhere 


Af(kT)  =  (1  -  f)  Af{(k-1)T}  +  ?  A(kT) 


A(kT)  =  input(output)  difference  signal 
Af(kT)  =  filtered  input(output)  difference  signal 
?  =  filter  constant 


(4-15) 


The  filter  constants  must  be  identical  in  order  to  preserve  the 
input-output  relationship  that  defines  the  plant.  The  value  of  ?  is 
selected  again  by  trial  and  error  as  a  tradeoff  between  the  desire  to 
attenuate  noise  spikes  in  the  input-output  data,  and  the  need  to 
maintain  sufficient  excitation  in  the  input-output  perturbation  signals. 
The  value  of  =  0.2  represents  an  acceptable  compromise  for  this 
study. 

To  deal  directly  with  the  problem  of  noise  spikes,  and/or  abrupt 
variations  of  the  parameter  estimates  in  general,  the  estimates 
the~selv'es  are  filte*'ed  before  they  are  used  for  any  control  law 
calculation.  The  for'm  of  the  low  pass  filter  chosen  for  this  function 
is  a  discreticed  version  of  a  simple  first  order  analog  filter  given  by 


where  s  is  the  Laplace  operator,  T  is  the  sampling  period  and  z  is 
the  discrete  Z  transform  operator.  This  results  in  the  following 
filter  algorithm 


Qf(kT)  =  Q  Gf{(k-1)T}  +  .-2  [  G(kT)  +  G{(k-1)T}  ]  (4-18) 


with 


<■1  =  2  :  r 

^'2  =  2  T~  (^-20) 

where  Gf(kT)  is  the  filtered  parameter  estimate. 

Although  this  low-pass  filtering  helps,  often  the  transients  are  of 
such  magnitude  that  the  filter  may  "charge-up"  to  a  very  large  offset  as 
compared  to  the  magnitude  of  the  correct  parameter  estimate.  In  that 
case  it  may  take  the  filter  a  Icng  time  to  return  to  the  nominal 
estimate  depending  on  the  value  chosen  for  v,  even  though  the  estimator 
itself  may  return  quickly  to  the  correct  value.  This  situation  can 
cause  a  slower  rate  of  adaptation  than  desired.  A  very  simple  way  of 
solving  the  problem  is  to  implement  a  rate  limiter  at  the  input  of  the 
parameter  filter  to  restrict  the  magnitude  of  the  fluctuations  in  the 
incoming  estimates.  Because  the  relative  magnitude  of  ti.e  individual 
parameters  may  vary  significantly,  the  rate  limit  is  not  set  to  an 
absolute  quantity,  but  rather  expressed  as  a  percentage  variation  from 


the  last  input  received  to  the  parameter  estimate  filter.  The  final 
version  of  the  parameter  estimate  filter  is  then  given  by  the  following 
algorithm 

For  each  element  Q''  of  the  parameter  vector  e,  iterate  thru 


IF  QT{(k-l)T}  >  A  THEN 

(4-21a) 

A0iRL(kT)  =  Q-i{(k-l)T}| 

(4-21b) 

Q-iRL(kT)  =  GlRL{(k-l)T}  +  ... 

MAx|  -  A0iRL(kT)  ,  MIN[(0i(kT)  -  0iRL{ (k-1 )T} )  ,  AO^RLCkT)  ] 

(4-21C) 

ELSE 

0TRL(kT)  =  0i(kT) 

(4-21d) 

END  IF 

oV(kT)  =  QV{(k-l)T}  +  ^2  [  Q^RlCkT)  +  O^RL{(k-l)T}  ]  (4-21e) 


where 

A  is  a  very  small  number  (  =  le-6) 

AC'’RL(kT)  is  a  rate  limited  parameter  estimate  increment 

I  is  the  allowable  percentage  of  the  last  input  to  form  the 
rate  limit 

o^’RL(kT)  is  a  rate  limited  parameter  estimate 


and  e''f(kT)  ,  (i  ,  (2  defined  as  before.  The  niagnitude  test  in  the 
begining  IF  statement  is  necessary  in  case  the  parameter  estimate  passes 
through  zero.  If  this  test  were  not  included,  because  the  maximun  and 
minimun  limits  in  Eqn  (4-21c)  are  dependent  on  the  magnitude  of  the  last 
incoming  estimate,  the  filter  would  "lock  on"  to  a  value  of  zero. 
Therefore,  when  the  estimates  are  within  a  very  small  band  around  zero 
the  limiter  is  bypassed  and  the  estimates  are  fed  directly  into  the 
low-pass  filter  as  Eqn  (4-21d)  indicates. 

The  filter  pole  location  and  the  input  rate  limiter  percentage 
value  z  are  both  design  parameters  and  are  adjusted  by  trial  and 
error.  The  values  of  o.-  =  2.25  and  i  =  25  are  used  throughout  the 
study.  The  filtering  procedure  just  discussed  helps  in  preventing 
transients  to  be  passed  on  to  the  calculation  of  the  control  law  in  Eqns 
(2-47)  and  (2-48)  and  therefore  improves  the  performance  of  the  adaptive 
control  system. 

Another  way  of  improving  the  the  performance  of  the  estimation 
algorithm  in  the  presence  of  noise  is  by  proper  scaling  of  the  parameter 
vector  0  and  measurement  matrix  v.  Although  not  formally  a  signal 
filtering  technique,  the  scaling  of  these  variables  helps  to  increase 
the  dynamic  range  of  the  identification  procedure  in  regards  to  noise 
level  effects.  When  the  standard  deviation  of  the  noise  level  is 
roughly  the  same  order  of  magnitude  of  some  of  the  elements  in  the 
parameter  vector,  the  estimates  of  those  parameters  may  suffer  from 
large  transients  and  biases  that  may  lead  to  close-loop  instability. 
This  is  specially  true  for  paramete''s  that  are  very  close  to  zero,  for 
example,  elements  of  the  step-response  matrix  corresponding  to  a  model 


of  an  aircraft  at  a  flight  condition  characterized  by  low  dynamic 
pressure  and  low  control  surface  effectiveness.  If  an  estimate  of  a 
particular  element  of  H(T)  goes  to  zero  or  oscilates  substantially 
around  that  region  because  of  noise  or  poor  input-output  excitation,  the 
control  law  gains  become  undefined  or  vary  wildly  because  of  the  inverse 
relationship  between  H(T)  and  the  control  law  gains.  This  situation  may 
then  lead  to  instability  if  proper  precations  are  not  taken. 

The  scaling  of  G  and  T  may  be  used  to  help  reduce  the  effects  of 
transients,  and  it  can  be  carried  out  in  various  ways.  The  method  used 
here  comes  from  Eqns  (3-9)  and  (3-11).  If  one  considers  the  response  of 
the  dicrete  system  given  by  Eq  (3-9)  at  t  =  T  .with  zero  initial 
conditions  and  input  u(0),  we  have 


y(T)  =  Bi  u(0)  +  .(T) 


~  TCB  u(0)  +  t (T) 


(4-22a) 


(4-22b) 


Eqn  (4-22)  can  also  be  expressed  as 


y(T)  .  I  TCB  T  (y(0)  ,  +  -  (T) 


(4-23) 


which  suggest  that  by  multiplying  the  elements  in  o  by  the  sampling 
frequency  (f  =  scale  factor),  and  the  measurement  matrix  v  by  the 
sampling  time,  the  magnitude  of  the  parameters  can  be  scaled  away  from 
the  noise  le.el  while  maintaining  the  same  inpot-cutput  relationship. 
Other  scale  factors  may  be  applied,  hc^ever  ccnside^'ation  must  be  given 


always  to  avoid  making  the  magnitude  of  the  elements  of  measurement 
matrix  too  small.  The  scaling  procedure  makes  it  possible  for  the 
algorithm  to  provide  better  performance  in  environments  with  increased 
sensor  noise  by  reducing  the  effect  of  parameter  estimate  transients  in 
the  closed-loop  response  of  the  system. 

The  case  in  which  the  estimates  go  to,  or  pass  through  zero 
momentarily  are  handled  by  evaluating  the  step-response  matrix  before  it 
is  used  in  the  control  law  calculations.  Parameter  estimates  resulting 
in  a  singular  step-response  matrix  are  discarded  and  the  control  law 
gains  are  not  updated  until  new  estimates  are  obtained  which  result  in 
an  invertible  H(T). 

As  it  was  mentioned  in  section  4.3.2,  all  of  the  elements  composing 
the  identification  algorithm  and  control  law  design  equations  are 
implemented  in  MATRIX^  in  the  block  labeled  "ADAPT2"  (Figure  4-10).  A 
more  detailed  look  of  which  is  shown  in  Figure  4-70.  The  delay  blocks 
labeled  "DEL  U"  and  "DEL  Y"  generate  delayed  samples  of  the  input  and 
output  signals  used  in  forming  the  measurement  matrix  T.  The  block 
labeled  "ADAPT"  is  a  FORTRAN  subroutine  that  implements  the  estimation, 
filtering  algorithms  and  safety  checks  just  described,  and  the  block 
labeled  "P  INDX"  implements  the  calculation  of  the  model-following 
criteria  given  by  Eqn  (4-10)  as  shewn  in  Figure  4-71. 

To  consolidate  the  discussion  on  the  design  of  the  recursive 
identification  algorithm,  a  list  summarizing  the  values  of  the  design 
parematers  used  for  the  estimation,  and  the  signal  conditioning  filters, 
is  provided  in  Table  4-3.  Chapter  5  will  sum.marize  the  results  obtained 
comparing  the  performance  of  both,  fixed  gain,  and  adaptive  controllers 


Table  4-3 


Final  Estimation  Algorithm  Design  Parameters 


Parameter  Estimation  Algorithm 
a  =  5.0E-5 
Fault  Detector 

7]^  =  0.85 

72  =  0.95 
rg  =  0.5 

Noise  variance  estimator 

73  =0.95 

r  =  20 

ri  =  0.2 


Input/Output  Perturbation  Filter 

f  =  0.2 

Parameter  Estimate  Filter 

=  2.25 
z  =  25  % 


under  conditions  of  plant  parameter  changes  and  sensor  noise. 


5.1 


Introduction 


This  chapter  presents  the  simulation  responses  obtained  by 
applying  the  design  procedures  described  in  chapter  4  of  this  thesis  to 
an  AFTI/F-16  aircraft  model  used  as  representative  of  the  VISTA/F-16 
dynamics. 

The  main  objective  of  this  study  is  to  determine  the  effectiveness 
of  a  parameter-adaptive  control  system  in  maintaining  model-following 
fidelity  for  an  in-flight  simulator  under  changing  conditions.  This  is 
accomplished  by  properly  adjusting  the  control  law  based  on  information 
derived  from  input-output  data  measurements  that  reflect  the  current 
characteristics  of  the  host  vehicle. 

This  chapter  is  divided  into  two  main  sections  corresponding  to  the 
evaluations  performed  on  a  fixed  gain,  and  a  parameter-adaptive  control 
system.  Each  section  in  turn  includes  responses  detailing  the 
individual  system's  behavior  under  conditions  of  plant  parameter  changes 
and  sensor  noise.  This  is  done  to  establish  a  comparative  base  upon 
which  the  virtues  and  deficiencies  of  each  approach  can  be  pointed  out. 

5 . 2  Fived  Gain  Controller  Responses 

5.2.1  Plant  Parameter  chance.  To  assess  the  effect  of  either 
plant  parameter  uncertainties  or  changes  in  the  model  following  tracking 
performance,  a  fixed  gain  control  law  based  on  the  design  parameters  of 
Table  4-2  is  tested.  The  input  commands  of  Figures  4-12  and  4-13  are 


used.  At  six  seconds  into  the  simulation  the  plant  dynamics  are  changed 
from  a  model  at  Mach  0.9  to  one  at  Mach  0.3  (both  at  10,000  ft.),  to 
obtain  information  on  the  robustness  properties  of  the  control  law.  The 
simulation  responses  are  shown  in  Figures  5-1  thru  5-8. 

Figure  5-1  shows  both,  the  flight  path  command  input  (solid  line) 
superimposed  on  the  host  aircraft  response  (dashed  line).  It  is  obvious 
from  this  plot  that  the  control  law  offers  sufficient  robustness  to 
maintain  the  level  of  tracking  for  this  command  input.  This  is  also 
exemplified  in  Figure  5-2,  which  shows  the  model-following  criteria  as 
calculated  from  Eqn  (4-10).  The  solid  line  represents  the  average  error 
absolute  value  response,  and  the  dashed  line  represents  the  average 
absolute  value  response  of  ten  percent  of  the  reference  signal.  The 
error  response  is  clearly  below  the  performance  limit  specification. 

In  the  case  of  the  pitch  rate  response,  ttie  situation  changes 
drastically  as  shown  in  Figures  5-3  and  5-4.  These  Figures  follow  the 
same  notation  as  Figures  5-1  and  5-2,  as  well  as  all  the  remaining 
aircraft  and  performance  index  plots  in  this  thesis.  Although  the 
response  of  Figure  5-3  demonstrates  the  ability  of  this  particular 
control  law  design  technique  to  provide  satisfactory  control  in  the 
presence  of  plant  parameter  uncertainty,  the  magnitude  of  the  change  in 
plant  dynamics  is  sufficient  to  violate  tne  model-following  performance 
criterion,  as  depicted  in  Figure  5-4.  Figures  5-5  and  5-6  show  the  host 
aircraft's  elevator  deflection,  and  elevator  deflection  rate.  Those 
corresponding  to  the  flaperon  deflection  and  deflection  rate  are  shown 
in  Figures  5-7  and  5-3.  All  of  these  control  surface  responses  are  well 
within  bounds,  and  considered  acceptable. 
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Figure  5-1.  Flight  path  angle  command  and  response  (deg) 
Fixed  gain  control  law.  Plant  parameter  change. 
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Figure  5-3.  Pitch  rate  command  and  response  (deg/sec) 
Fixed  gain  control  law.  Plant  parameter  change. 


Figure  5-4.  Pitch  rate  tracking  performance  index  (deg/ 
Fixed  gain  control  law.  Plant  parameter  change. 


Figure  5-5.  Elevator  deflection  (deg). 
Fixed  gain  control  law.  Plant  parameter  change 


Fig'jre  5-7.  Flaperon  deflection  (deg). 
Fixed  gain  control  law.  Plant  parameter  change 


Figure  5-S.  Flaperon  deflection  rate  (deg, 'sec) 
Fixed  gain  control  law.  Plant  parameter  change 


5.2.2  Plant  Parameter  Change  and  Sensor  Noise.  In  addition  to 
robustness  to  plant  parameter  changes,  it  is  important  to  demonstrate 
the  influence  of  sensor  measurement  noise  in  the  closed-loop  system 
response.  Each  plant  output  measurement  is  corrupted  by  zero-mean, 
white  gaussian  noise  added  to  the  simulation  by  the  "noise"  block 
location  in  Figure  4-8.  Figures  5-9  thru  5-16,  and  5-17  thru  5-24  show 
two  representative  samples  of  the  fixed  gain  system  responses  to  noise 
levels  with  standard  deviations  of  0.0573  deg  (deg/sec  for  pitch  rate), 
and  0.314  deg  (deg/sec  for  pitch  rate)  respectively. 

The  most  noticeable  effect  of  sensor  measurement  noise  is  that  of 
increased  control  surface  deflection  and  deflection  rate  activity, 
particularly  the  latter.  The  effect  on  tracking  performance,  however, 
is  hardly  noticeable. 

Although  the  control  surface  rates  shown  in  these  Figures  are  below 
the  specified  limits,  the  sustained  level  of  effort  required  from  the 
actuators  indicated  by  these  Figures  may  represent  a  severe  demand  on 
the  aircraft's  hydraulic  system.  These  results  merely  reiterate  the 
need  for  high  quality  sensors,  in  combination  perhaps,  with 
complementary  filter  functions  to  achieve  low  noise  level 
specifications.  This  is  particularly  appropriate  in  cases  where  high 
ga'n  control  laws  like  the  ones  usee  in  this  thesis  are  applied.  The 
need  for  lew  sensor  noise  levels  will  become  mere  evident  when  the 
inp^t-output  measurements  are  used  to  update  the  control  law  gains  as  in 
the  case  of  the  parameter  adaptive  control  law  which  is  presented  next. 


Figure  5-9.  Flight  path  angle  command  and  response  (deg). 
Fixed  gain  control  law.  Plant  parameter  change  and  sensor  noise 

(std.  dev.  =  0.0573  deg) 


FigL.'re  5-10.  Flight  path  angle  tracking  performance  index  (deg) 
Fixed  gain  control  law.  Plant  parameter  change  and  sensor  noise 

(std.  dev.  =  0.0573  deg) 


5.3 


5.3.1  Fixed  Plant  Dynamics.  To  becoir,e  familiar  with  the 
characteristics  of  the  parameter-adaptive  system  under  various 
conditions,  its  performance  is  tested  in  several  stages.  The  adaptation 
mechanism  is  first  tested  with  simulations  in  which  the  dynamics  of  the 
aircraft  model  (at  MACH  =  0.9  ,  10,000  ft.)  are  kept  constant.  The 
purpose  of  this  test  is  simply  to  verify  the  capability  of  the 
identification  algorithm  to  track  the  parameters  of  a  known  model. 
Initial  testing  is  done  without  scaling  the  difference  equation  model 
(see  section  4.4.5).  The  identification  algorithm  is  initialized  in  all 
simulations  with  the  parameter  vector  obtained  from  the  transfer  matrix 
of  the  indicated  model  (see  appendix  A).  The  covariance  matrix  is  set 
to  a* I,  and  the  prediction  error  variance  is  set  to  1  x  lO'^O  deg^ 
(deg2/sec2).  The  adaptation  mechanism  is  not  switched  on  until  two 
seconds  into  the  simulation  to  allow  time  for  the  measurement  matrix  to 
fill  up  with  past  input-output  data  and  for  the  prediction  error 
variance  calculation  (Eqn  3-60)  to  stabilize. 

Figures  5-25  thru  5-32  show  the  aircraft  responses  in  this 
simulation.  As  expected,  these  do  not  show  variation  in  tracking 
performance  from  those  of  the  nominal  design  case  of  Figures  4-22  thru 
4-29.  The  only  noticeable  difference  is  a  slight  increase  in  control 
surface  deflection  rates.  Tnis  is  attributed  to  the  transients  that 
occur  in  the  filtered  parameter  estimates  that  are  used  to  update  the 
control  law  gains.  Figure  5-33  shows  the  fault  detector  test  signal 
r(kT).  Although  the  signal  snows  several  significant  peaks,  it  remains 
at  all  times  below  the  selected  threshold  of  cq  =  0.5,  thus  indicating 
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Figure  5-25.  Flight  path  angle  coiruiand  and  response  (deg). 
Adaptive  control  law.  No  plant  parameter  change. 
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Figure  5-26.  Flight  path  angle  tracking  performance  index  (deg). 
Adaptive  control  law.  No  plant  parameter  change. 
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Figure  5-27.  Pitch  rate  command  and  response  (deg/sec). 
Adaptive  control  law.  No  plant  parameter  change. 


gure  5-28.  Pitch  rate  tracking  performance  index  (deg/sec). 
Adaptive  control  law.  No  plant  parameter  change. 
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Figure  5-29.  Elevator  deflection  (deg). 
Adaptive  control  law.  No  plant  parameter  change 
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Figure  5-30.  Elevator  deflection  rate  (deg/sec) 
Adaptive  control  law.  No  plant  pat'ameter  change 


i 


Figure  5-31.  Flaperon  deflection  (deg). 
Adaptive  control  law.  No  plant  parameter  change 


Figure  5-32.  Flaperon  deflection  rate  (deg/sec) 
Adaptive  control  law.  No  plant  parameter  change 
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Figure  5-33.  Es 


mate  of  step-response  matrix  element  h2i(kT) 


Figure  5-39.  Covariance  matrix  element  P33(kT). 


that  no  abrupt  plant  parameter  changes  are  reported.  This  is  also 
expexted  since  no  change  in  aircraft  dynamics  has  taken  place.  Figures 
5-34  thru  5-41  show  the  estimates  of  the  elements  of  the  step-response 
matrix,  both  filtered  (solid  lines)  and  unfiltered  (dashed  lines),  and 
the  cori'espcnding  diagonal  elements  of  the  parameter  covariance  matrix. 
Although  these  plots  show  the  estimates  undergoing  several  transients, 
their  average  values  remain  close  to  the  actual  parameters.  The 
variations  in  the  magnitude  of  the  filtered  estimates  are  considered 
tolerable  given  the  robustness  properties  of  the  control  law. 

The  periods  in  which  the  estimates  show  the  greater  deviations  are 
usually  associated  with  the  instances  where  the  control  surface 
deflection  rates  undergo  spikes.  These  spikes  are  due  to  abrupt  changes 
in  magnitude  and/or  direction  in  the  command  signals.  Here  in 
particular  those  of  the  pitch  rate  com.mand.  When  the  estimation 
algorithm  detects  the  higher  level  of  excitation  to  the  plant,  it 
quickly  discounts  old  data  to  keep  up  with  possible  changes  in  the  plant 
dynamics.  However,  the  sudden  removal  of  that  excitation  (surface 
deflection  rate  spike  ends)  causes  the  algorithm  to  slow  down  (or  step) 
momentarily  its  discounting  of  old  information  to  match  whatever 
infermatien  content  the  current  input-output  data  might  bring.  This 
also  repT'ese"ts  a  slew  down  in  adaptation  rate.  This  action  may  happen 
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to  minimize  the  least-squares  cost  function. 


That  being  f  e  case,  it  ma_.-  take  seme  time  for  tne  algorithm  to  acquire 
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simulation  setup  is  used  here  again,  this  time  to  assess  the  effects  of 
sensor  measurement  noise  in  the  identification  pr-ccedure.  Because  of 
the  relatively  small  magnitude  of  some  of  the  elements  of  the 
step-response  matrix,  the  noise  level  that  can  be  tolerated  (without  any 
additional  safety  net  type  features)  is  significantly  smaller  than  that 
of  the  fixed  gain  system.  Figures  5-42  thru  5-53  show  the  results 
obtained  with  a  noise  standard  deviation  of  0.00181  deg  (deg/sec). 

The  aircraft's  responses  in  this  case  are  almost  identical  to  those 
in  the  previous  simulation  with  the  adaptation  mechanism,  but  with  still 
higher  control  surface  deflection  rates  due  to  the  corrupted 
measurements  and  transients  in  the  parameter  estimates  discussed 
earlier.  Despite  this,  no  degradation  of  tracking  performance  occurs. 
In  contrast  to  the  tracking  performance,  the  responses  of  the 
identification  algorithm  are  significantly  different  from  the  previous 
ones.  The  fault  detector  signal,  for  example,  is  shifted  in  mean  value 
from  zero  to  -0.5.  The  reason  for  this  behavior  is  that  the  fault 
detector  algorithm  is  derived  under  the  assumption  of  an  equation  error 
term  with  "white"  noise  characteristics  around  the  nominal  parameter 
estimates.  The  introduction  of  the  input-output  perturbation  filter 
given  by  Eqn  (4-15)  has  the  side  effect  of  reducing  the  noise  spectrum, 
thus  producing  data  corrupted  bg  "colored"  noise  instead.  It  is  a  well 
known  fact  that  "colored"  noise  produces  biases  in  the  parameter 
estimates,  and  it  is  those  biases  which  ai'e  being  reflected  in  the  mean 
of  the  fault  detector  signal.  It  is  impoi'tant  to  note  however,  that  the 
property  of  the  test  signal  as  an  indicator  of  abrm.rt  plant  param.eter 
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Figure  5-42.  Flight  path  angle  command  and  response  (deg). 
Adaptive  control  law.  No  plant  parameter  change  and  sensor  noise. 

(std.  dev.  =  0.00181  deg). 
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Figure  5-43.  Flight  path  angle  tracking  perfcrmar^ce  ir,:-.- 
Adaptive  control  law.  fio  plant  parameter  chanoe  a-  1  s-r*  ■ 

(std.  dev.  =  O.CGISI  deg).’' 
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Figure  5-44.  Pitch  rate  command  and  response  (deg/s 
Adaptive  control  law.  No  plant  parameter  change  and  sens 


5-45.  Pitch  rate  tracking  performance  index  (deg/sec), 
control  law.  No  plant  parameter  change  and  ser^':"r  noise, 
(std.  dev.  =  0.00181  deg/sec). 
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Figure  5-46.  Elevator  deflection  (deg). 

Adaptive  control  law.  No  plant  parameter  change  and  sensor  noise 
(std.  dev.  =  0.00181  deg  (deg/sec)). 
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Figure  5-47.  Elevator  deflection  rate  (deg/sec). 
Adaptive  control  law.  No  plant  parameter  change  and  sensor  noise 
(std.  dev.  =  0.00181  deg  (deg/sec)) 
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Figure  5-48.  Flaperon  deflection  (deg). 

Adaptive  control  law.  No  plant  parameter  change  and  sensor  noise 
(std.  dev.  =  0.00181  deg  (deg/sec)). 


Figure  5-49.  Flaperon  deflection  rate  (deg/sec). 
Adaptive  control  law.  No  plant  parameter  change  and  sensor  noise 
(std.  dev.  =  O.OOISI  deg  (deg/sec)). 
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Figure  5  50.  Fault  detector  test  signal  with  sensor  noise, 
(std.  dev.  =  0.00181  deg  (deg/sec)). 
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Figure  5-55.  Estimate  of  step-response  matrix  element  hp^kT)  with 
sensor  noise,  (std.  dev.  =  0.00181  deg  (deg/sec)). 
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Figure  5-56.  Covariance  matrix  element  pggCkT)  with  sensor  noise, 
(std.  dev,  =  O.OOISI  deg  ('deg/sec)). 
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Figure  5-57.  Estimate  of  step-response  matrix  element  h22(kT)  with 
sensor  noise,  (std.  dev.  =  0.00181  deg  (deg/sec)). 
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Figure  5-53.  Covariance  matrix  element  P44(kT)  with  sensor  noise, 
(std.  dev.  =  0.00181  deg  (deg/sec)). 
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changes  remains  unchanged.  The  shift  in  the  mean  value  of  r(kT)  is 
easily  compensated  for  by  simply  adjusting  the  threshold  levels  of  ro 
and  ri  to  account  for  the  new  mean  signal  value  of  -0.5.  This  action 
takes  place  a  few  samples  after  the  algorithm  is  started,  to  allow  for 
initalization  transients  to  fade  out. 

The  parameter  estimates  in  this  simulation  show  considerable  noise 
jitter  superimposed  onto  them  with  a  corresponding  increase  in  their 
uncertainty  as  signified  by  the  plots  of  their  respective  covariance 
matrix  diagonal  elements.  The  fluctuations  of  these  estimates  around 
their  nominal  values  however,  are  significantly  smaller  than  those  of 
the  simulation  without  noise.  The  reason  for  this  is  the  additional 
excitation  that  the  plant  is  receiving  from  the  control  surfaces,  which 
helps  the  identification  process.  This  additional  control  surface 
activity  is  the  result  of  the  control  law  trying  to  maintain  tracking 
performance  despite  the  fluctuations  of  the  measured  variables  due  to 
sensor  noise. 

5.3.3  Plant  Parameter  Change.  The  next  step  in  the  testing  of 
the  parameter-adaptive  system  is  to  assess  its  capability  to  track  the 
changes  in  plant  parameters  and  consequently  maintain  the  tracking 
performance  at  a  specified  level.  This  is  accomplished  by  running  the 
simulation  and  specifying  a  change  in  the  plant  dynamics  model  to  occur 
at  six  seconds  into  the  simulation.  The  new  set  of  dynamics  corresponds 
to  those  of  the  AFTI/F-16  at  a  flight  condition  of  mach  0.9  at  10,000 
ft.  MSL  (Appendix  A).  The  simulation  responses  are  shown  in  Figures 
5-59  thru  5-75.  When  identification  of  the  step-response  matrix 
elements  takes  place,  the  tracking  performance  improves  dramatically  as 
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Figure  5-59. _  Flight  path  angle  command  and  response  (deg) 
Adaptive  control  law.  Plant  parameter  change. 
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Figure  5-60.  Flight  path  angle  tracking  performance  index  (deg) 
Adaptive  control  law.  Plant  parameter  change. 
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Figure  5-62.  Pitch  rate  tracking  performance  ind 
Adaptive  control  law.  Plant  parameter  ch 
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Figure  5-65.  Flaperon  deflection  (deg). 
Adaptive  control  law.  Plant  parameter  change. 
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Figure  5-65.  Flaperon  deflection  rate  (deg/sec). 
Adaptive  control  law.  Plant  parameter  change. 


.'•.V 


■'V/ 


V'  j' 


B  10  ir 


T-.r-E  (  SEC  ) 


IB  13 


Figure  5-68.  Estimate  of  step-response  matrix  element  hii(kT). 

Plant  parameter  change. 
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Figure  5-69.  Covariance  matrix  element  p]^i(kT), 
Plant  param.eter  change. 


Figure  5-70. 


Estimate  of  step-response  matrix  element  hi2(l<T) 
Plant  parameter  chnage. 


Figure  5-72.  Estimate  of  step-response  matrix  element  h2i(kT) 

Plant  parameter  chnage. 
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Figure  5-73.  C: . ar iance  matrix  element  P33(kT). 
Plant  para"'ete''  chanae. 


Figure  5-74.  Estimate  of  step-response  matrix  element  h22(kT) 

Plant  parameter  change. 
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Figure  5-75.  Covariance  matrix  element  P44(kT). 
Plant  parameter  change. 


shown  in  Figures  5-61  and  5-62  for  the  pitch  rate  response.  This  is  in 
comparison  to  the  performance  of  the  fixed  gain  system  in  Figures  5-3 
and  5-4. 

Another  noticeable  characteristic  in  these  aircraft's  reponses  is 
the  encountering  of  flaperon  deflection  rate  limits.  This  is  attributed 
to  the  change  in  direction  of  the  pitch  rate  command  at  those  particular 
instants,  in  combination  with  fluctuations  suffered  by  the  parameter 
estimates  at  the  new  flight  condition.  The  new  set  of  plant  dynamics  is 
characterized  by  a  step-response  matrix  with  significantly  smaller 
elements  than  the  previous  one.  Fluctuations  in  these  small  numbers  are 
bound  to  produce  large  variations  in  the  control  law  gain  matrices,  and 
therefore  larger  deflection  rates.  The  encountering  of  flaperon  rate 
limits  however,  is  not  considered  a  severe  problem  in  this  simulation 
since  the  instances  in  which  the  flaperon  is  rate  limited  are  very 
short,  therefore  not  affecting  signif icnatly  the  stability  of  the 
system. 

The  performance  of  the  identification  algorithm  is  considered  very 
good.  The  fault  detection  algorithm  in  this  simulation  clearly 
indicates  the  abrupt  change  in  the  plant  dynamics  as  shown  in  Figure 
5-67  by  the  peak  of  r(kT)  exceeding  the  threshold  rg  =  0.5.  The 
abrupt  change  in  the  plant's  step-response  matrix  elements,  together 
with  the  identification  of  the  new  parameters  and  parameter  variances  is 
shown  in  Figures  5-63  thru  5-75. 


5.3.4  Plant  Pa*~a'''eter  Chance  a^d  Se'^scr  Ncise.  To  present  a 
worst  case  scenario  to  the  adaptation  mecnanism,  the  si~jlati:n  with 


abrupt  plant  parameter  changes  is  r.Zf>  e-ec„ted  with  the  add'ticn  c' 


Figure  5-78.  Pitch  rate  coirmand  and  response  (deg/sec). 
Adaptive  control  law.  Plant  parameter  change  and  sensor  noise 
(std.  dev.  =  0.00181  deg/sec). 
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Figure  5-82.  Fiaperon  deflection  (deg). 

Adaptive  control  law.  Plant  parameter  change  and  sensor  noise 
(std.  dev.  =  0.00181  deg  (deg/sec)). 


Figure  5-S3.  Flacercn  deflection  rate  (deg/sec). 
Adaptive  control  law.  Plant  para"’eter  change  a",1  sensor  noise 
(std.  dev.  =  O.OCISI  deo  (deo'seol). 
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Figure  5-85.  Estimate  of  step-response  matrix  element  hii(kT)  with 
sensor  noise,  (std.  dev.  =  0.00181  deg  (deg/sec)). 
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Figure  5-85.  Covariance  matrix  ele-^ent  pi’(kT)  ;vith  sensor  noise, 
(std.  dev.  =  O.OOISI  deg  (oeg/sec)). 


Figure  5-87.  Estiriiate  of  step-response  matrix  element  hi2(kT)  with 
sensor  noise,  (std.  dev.  =  0.00181  deg  (deg/sec)). 


Figure  5-83.  Covariance  matrix  element  p-?(kT)  with  sensor'  noise, 
(std.  dev.  =  0.00181  deg  Cceg/sec)). 


sensor  measurement  noise.  The  system's  performance  in  this  simulation 
is  very  similar  to  that  of  the  previous  one  until  the  last  few  seconds, 
where  the  closed-loop  system  develops  an  instability  in  the  aircraft's 
responses.  The  instability  is  the  result  of  control  surface  rate 
limiting.  This  condition  is  caused  in  turn  by  noise  induced 
fluctuations  in  the  filtered  parameter  estimates,  at  a  time  when  these 
estimates  are  extremely  close  to  the  zero  axis.  As  it  was  mentioned 
earlier,  the  situation  just  described  causes  large  changes  in  the 
control  law  gains  which  are  responsible  for  driving  the  control  surfaces 
to  their  motion  rate  limits. 

To  alleviate  this  problem  the  simulation  is  then  executed  with 
scaled  parameter  vector  and  measurement  matrix  as  described  in  section 
4.4.5.  The  results  are  ilustrated  in  Figures  5-93  thru  5-109.  The 
responses  of  the  closed-loop  system  are  stable  and  tracking  performance 
is  maintained  as  desired.  The  performance  of  the  identification 
algorithm  is  clearly  better  when  scaled  variables  are  used.  The  fault 
detector  test  signal  in  Figure  5-101  shews  a  prominent  peak  shortly 
after  the  time  of  the  plant  dynamics  change,  indicating  unmistakably  the 
occurence  of  the  fault.  The  parameter  estimates  themselves  exhibit  a 
significantly  smoother  behavior  (Figures  5-102,  5-104,  5-106,  and 
5-103),  thus  the  filtered  estimates  are  virtually  free  of  oscillations 
or  fluctuations  around  the  zero  axis.  This  behavior  accounts  for  the 
improved  performance  of  the  clcsed-locp  system  responses  over  the  ones 
shewn  using  the  unsealed  algcrithm,  for  the  level  of  noise  used. 

The  scaled  version  of  the  algorithm  is  again  tested,  this  time  with 
a  higher  noise  level  of  standard  deviation  of  O.GG5'’3  deg  (deg/sec). 
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Figure  5-93.  Flight  path  angle  command  and  response  (deg). 
Adaptive  control  law.  Plant  parameter  change  and  sensor  noise, 
(std.  dev.  =  0.00181  deg).  Scaled  algorithm  (SF  =  100). 


Figure  5-94.  Flight  path  angle  tracking  performance  inCe<  (de 
Adaptive  control  law.  Plant  parameter'  change  and  senscr  nois 
(std.  dev.  =  O.OOISI  deg).  Scaled  algcrith""  (SF  =  ICO). 
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Figure  5-103.  Estimate  of  step-response  matrix  element  h22(t<T)  with 
sensor  noise,  (std.  dev.  =  0.00181  deg  (deg/sec)). 

Scaled  algorithm  (SF  =  ICO). 
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The  closed-loop  system  response  still  exhibits  a  satisfactory  behavior, 
although  with  noticeably  higher  levels  of  flapercn  deflection  rates 
(Figure  5-117).  The  bursting  type  of  behavior  shown  by  the  flaperon 
deflection  rate  between  seven  and  eight  seconds  in  the  simulation  has 
its  origin  in  the  fluctuation  of  the  filtered  estimate  of  h22('<T)  close 
to  the  zero  axis  during  the  same  period  of  time.  Overall,  the 
performance  of  the  identification  algorithm  is  good  with  considerably 
less  noise  superimposed  on  the  estimates,  although  the  level  of  noise  in 
this  simulation  seems  to  be  approaching  the  limit  for  which  adequate 
protection  (against  fluctuations  close  to  the  zero  axis)  can  be  provided 
with  the  scaling  used  (scale  factor  =  100  =  (1/T)). 

To  verify  if  any  additional  beneficts  can  be  obtained  with  increased 
scaling,  a  simulation  was  executed  with  a  scale  factor  of  125.  All  the 
responses  were  identical  from  those  of  the  simulation  done  with  the 
scale  factor  of  100.  All  of  the  simulations  executed  with  sensor  noise 
seem  to  suggest  that  additional  safeguards  are  needed  in  the  case  of 
step-response  matrices  with  elements  that  are  very  close  to  zero. 
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5.3.5  Identification  of  Entire  Parameter  Vector  with  Full,  and 
Reduced  Order  Models.  As  a  final  test  of  the  adaptation 
mechanism,  the  estimation  of  the  plant  dynamics  is  carried  cut  letting 
all  cf  the  elements  of  the  parameter  vector  to  "float"  i.e.  estimating 
all  the  parameters  in  the  vector  diference  equation  model  cf  Eqns  (3-9) 
and  (3-10).  This  was  dene  for  a  full  order  model  (  N  =  4  ),  and  reduced 


order  models  of  third  (  N  =  3)  and  second  (  N  =  2  )  order.  All  of  the 
tests  resulted  in  unstable  aii'craft  responses.  Ttiis  is  because  the 


parameter  estimates  undergo  larce  transients,  particularly 
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Figure  5-116.  Flaperon  deflection  (deg). 

Adaptive  control  law.  Plant  para~eter  change  and  sensor  noise, 
(std.  dev.  =  0.QQ573  deg  (deg/sec)).  Scaled  algorithm  (SF  =  100). 
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Figure  5-117.  Flaperon  deflection  rate  (deg/sec). 

Adaptive  control  law.  Plant  parameter  change  and  sensor  noise, 
(std.  dev.  =  O.CG573  deg  (deg/sec)).  Scaled  algorithm  (SF  =  ICO). 
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Figure  5-125.  Estirate  of  step-response  matrix  element  h22(*<T)  with 
se'sor  noise,  (std.  dev.  =  0.CC573  deg  (deg/sec)). 

Scaled  algorithm  (SF  =  100). 
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aro'^nd  zero.  As  eAp.ainea  before,  these  transients  drive  the  control 
surfaces  to  their  deflection  rate  limits  for  too  long.  This  results  in 
an  eAcesive  phase  lag  for  the  control  action,  and  the  encountering  cf 
surface  deflections  limits.  All  of  these  factors  s..tsequently  result  in 
loss  of  conti'ol  of  tne  cpen-lcop  unstable  plant. 

The  large  transients  and  the  seemingly  long  convergence  time  for  the 
estimates  are  mainly  the  result  of  t^o  conditions.  First,  the  larger 
number  of  parameters  to  be  identified  puts  additional  demands  cn  the 
level  of  e/citation  reguired  fcr  gcod  identifcaticn  to  take  place  (see 
section  3.3).  With  an  evcitaticn  level  smaller  than  the  ideal  (for  the 
particular  r.-ber  cf  para~eters  involved)  the  estimator  will  tar.e  a 
longer  time  to  conv6'-ge  sc  tnat  it  can  accu-^.ulate  encuch  infcr-aticn 
that  allc.vs  the  calc..laticn  cf  a  reascnadj  precise  estimate. 

Final  Ij,  tne  pa- y-'vtricaticn  osed  f:-  defin-^g  the  inp..t-CMtput 
’’^oiel  s..ffe-s  from  s:  ca"eo  idcnti' iao  i  1  i  t^  pr. clems  en-e-g'ng  frc" 
pole-cer'o  ca-;elatic"  e'fects.  Wne-^  cotaining  the  :..tp-t  fL-nctions  for 
the  m^lti-inp:.t  mul t ■ -c..tp-t  plant  using  t:'e  c:“m:n  de'cminatcr  cbtained 
frc~  the  tra-'er  matri.,  fe  n;.~erat:r  pcl^ncmia's  ma^  be  fenced  tc 
contain  additional  te-ms  in  order  to  ratoh  the  fcr~  of  the  de^omina*-'!' . 
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Figure  5-142.  Estimate  of  step-response  matrix  element  h22(kT) 
Identification  of  full  4th  order  model. 


Figure  5-143.  Covariance  matrix 
Identification  of  full  4“h 


element  P44(I<T). 
order  model . 
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Figure  5-143.  Elevator  deflection  (deg). 
Adaptive  control  law.  Plant  parameter  change. 
Identification  of  full  order  model.'" 


Figure  5-149.  Elevator  ceflecticn  rate  (oeg  sec), 
Adaptive  control  law.  Plant  parameter  cha'’,ge. 
Identif ication  of  full  3’'4  cr-oe’'  r:c-l. 
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e  5-163.  Pitch  rate  corrin^.and  and  response  (deg/sec) 
Acapti'.e  control  law.  Plant  parameter  change. 
Identification  of  full  2^^  order  model. 


5-16d.  Fitch  rate  tracking  pct'fcr'mance  incex 
Adaptive  control  law.  Plant  panai^eter'  chant 
I'dentif icaticn  of  full  2"-"  cruer  mjcel. 
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jr'e  5-170.  Estimate  of 
Identif  icat^an 


step-respor'.s 
of  fall  2"''- 


e  natrix  element  hii(kT) 
cr'der  riicdel . 


.re  5-171.  Ca- a- 
Idertir icaticn 


full  2’'- 


contains  more  parameters  than  the  minirrio 


necessary  to  descrite  it, 


thet'efore  it  requires  ccnsideralole  more  excitation  for  proper' 
idtrotif icaticn.  If  identif icaticn  cf  all  tne  parao.ete’'s  in  the  rodel  cf 
toe  plant  is  necessary,  a  uniquel,.-  identifiacle  ridel  fcrr  s^cn  as  a 
ca-'iniial  autor'eor'ess  i  ve  roving  average  (A-t'A)  differ'e^oe  eo-aticr' 
S'':^lo  be  used  (11,  17,  Ic).  Given  that  tne  criginal  scope  cf  tnis 
ef'crt  was  to  identify  only  the  ele'"er,ts  of  the  matrix  coefficient  E]^, 
tne  selection  of  an  ARMA  form  or  a  transfer  matrix  derived  model 


le.els 


the 


its  tolerance  to  higher  noise  le.els,  the  reo..ction  of  para~eter 
estimate  transients  anO  their  effects  on  clcseh-lccp  perfcrrance. 


c.2.1  C^-^rol  Law  f':  iif icati-n.  Tne  si"-jlaticn  respcnses  presents 


in  AppenOix  D  shew  that  significant  irp--o.'e"-er,t  can  te  achieved  by 
i-p ting  anti-wind^p  previsions  in  t‘'e  c:''t’'ol  la.-..  Tne  particalan 
technique  used  in  this  in.-est i gat icn  is  based  si~plj  on  stepping  the 
calculaticn  of  the  integral  of  tne  e'Tcr'  in  the  cc^trol  law  whe'',e.er 
control  sur'fdce  (deflecticn  lirits  a’'e  r-a:hed.  Hc.-.'r.rr,  additicnal 
beneficts  may  be  obtained  by  modifying  tne  ccntr'ol 'er '  s  integral  action 
in  a  different  manne'".  Cons io-c-r' i ng  that  E:".s  (Z-:6)  and  (2-46)  e-press 
t''e  control  law  as  a  state-space  system  witn  a  pla"t  mat-ix  A  -  I 
(i.e.  w'th  cpen-lcop  poles  at  t‘'v  u-it  c'ncle),  it  St—'os  reascnatle  to 
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Most  of  the  problems  noted  in  the  simulations  have  their  root  in  the 
snail  magnitude  of  the  elements  of  the  step-respcne  natrix  (H(T)),  anp 
the  in.erse  relationship  between  H(T)  and  the  control  la.v  gain 
natr'ces.  For  that  reason,  it  is  ir^per'af  i . e  to  axoid  f  IcctL-aticn  the 
para-eter  estinates  v-hen  updating  tne  control  la.v  gains  since  t"is  ra. 


lead  to  higher  cortrol  s..!'fc 
instaPi lities. 
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Possible  corrective  neasL^-'es  that  can  te  tar.en  inon^de  the  additicn 
of  supervisory  type  functions  that  r.;nit:r  trer-ps  in  the  parareter 
esti~ates,  control  surface  de"  lectiC'S ,  anj  cjtp..'t  responses.  Updating 
of  t’'.e  control  la.v  gams  can  te  disc:-: •''"...ed  if  tr'ansients  in  the 
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This  input-output  model  structure  involves  the  esti^'iticn  of  Ejl  v-hich 
contains  "large"  numters  as  compared  with  those  in  E^,  thus  premising 
im, proved  noise  level  tolerance,  and  less  flu:t„aticr3  in  t’:'.e  para":eter 
estimates  for  a  given  noise  level.  If  necessa''.,' ,  scaling  of  tne 


para'".etet's  and  measurements  may  still  te  G;ne  in  a  simi'a''  man-'er  to 
that  described  in  chapter  4. 

The  identif  iabi  1  ity  problems,  long  convergence  time,  a-'d  la-'ge 
fluctuations  in  the  parameter  estimates  encountered  v.hen  identification 


of  the  entire  parameter  vevtor  was  considered  can  be  reduced  by 
employing  a  canonical  ARMA  model  representation  for  the  plant  djnamics 
as  mentioned  in  chapter  5.  This  model  structure  provides  fer  the 
minimal  r.^"'de''  of  para~eters  required  to  rept'ese":  the  p'a-.t,  tn..s 
el  i-'-’cti'-g  t:-e  pole-ce’'j  ca-celaticn  e'f-rsts  tnat  cre^r's  i-'  t-e  vector- 
d  i  f :  c--:"  :e  ec-ati:n  -ccvl  d--i.-c  fro"'  t’'e  t:'a-c“e-  -cf'--.  “liO.  s'-:e 
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This  increase  in  estimator  gain  may  be  appropriate  for  some  param.eters, 
but  excessive  for  others.  In  cases  when  a  large  number  of  parameters 
are  being  identified  but  only  a  few  may  be  changing  significantly,  this 
action  may  indeed  slow  down  the  identification  procedure  or  cause 
unnecessary  transients.  It  therefore  seems  reasonable  to  use  a  "vector" 
fault  detector  scheme  that  can  identify  changes  in  the  individual 
elements  of  the  parameter  vector.  Then,  the  appropriate  elements  of  the 
P  matrix  can  be  adjusted  in  a  manner  proportional  to  the  magnitude  of 
the  change. 

To  deal  with  the  problem  of  measurement  noise  directly,  the 
identification  algorithm  can  be  modified  to  employ  instrumental  varibles 
(IV)  and  extended  least  squares  (ELS)  techniques  (28)  that  are  better 
suited  to  handle  noisy  data  with  characteristics  different  from  those  of 
white  Gaussian  noise.  These  techniques  include  provisions  that  account 
for  the  shape  of  the  noise  spectrum,  or  correlate  out  the  effect  of 
noise  in  the  parameter  estimates  thus  reducing  biases  that  occur  when 
the  RLS  algorithm  is  used  in  a  "colored"  noise  environment. 

Further  extensions  to  this  research  should  consider  the 
modifications  named  above,  in  addition  to  study  problems  such  as  the 
identification  of  non-linear  models,  command  limiting  schemes,  control 
reconf iguraticn  strategies,  and  the  use  of  multiple-model  moving  bank 
estimator  techniques  (21).  It  was  assumed  for  this  investigation  that 
the  parameters  for  the  fixed  portion  of  the  difference  equation  model 
(representative  of  the  current  flight  condition)  were  available  in  an 
effort  to  reduce  the  number  of  parameters  that  needed  to  be  identified. 

.  atl'iihing  probabilistic  weightings  to  several  models,  it  may  be 


possible  to  select  the  best  "fixed  portion"  of  the  parameter  vector  for 
a  given  flight  condition.  Research  investigating  the  refinement  of  this 


technique  is  in  progress.  The  use  of  reduced  order  models  of  a 
canonical  ARMA  structure  also  needs  to  be  addressed  as  an  alternative  to 
multiple  models  scheme  proposed  above. 
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Appendix  A 

Aircraft  Data  for  State  Space  Models 


Tables  A-1  through  A-3  give  the  flight  pararieters,  aerodynamic 
derivatives,  and  diference  equation  models  for  the  flight  conditions 
used  in  this  thesis.  The  aerodynamic  data  was  obtained  from  a 
simulation  data  package  program  for  the  AFTI  F-16  aircraft,  available  at 
the  Flight  Dynamics  Laboratory's  Control  Synthesis  Branch.  The 
difference  equation  models  were  obtained  by  taking  the  aircraft  state 
space  models  and  generating  the  tranfer  matrix  from  their  discrete  time 
equivalent  representations.  This  was  done  with  the  aid  of  MATRIX^,  and 
its  model  manipulation  functions  "DISCRETIZE"  and  "TFORM"  (23). 


TABLE  A-1 

AFTI/F-16  Aircraft  Data 


Aircraft  Parameters 

c  (Wing  Mean  Aerodynamic  Chord)  =  11.32  ft. 
S  (Wing  Surface  Area)  =  300  ft^ 


b 

(Wing  Span) 

=  30  ft 

w 

(Weight)  =  21 

,018  lbs 

Moments  of  Inertia 

^xx 

=  14,145.68 

slugs-ft2 

^yy 

=  59,596.68 

slugs-f t^ 

I,, 

=  70,887.00 

slugs-ft2 

Ixz 

=  720.00 

slugs-ft2 
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Table  A-2 


Aircraft  Data  for  0.9  Mach,  10,000  ft 

IT  (dynamic  pressure)  =  825.33  Ib/in^ 

Vj  (trim  velocity)  =  969.66  ft/sec 
aj  (trim  angel  of  attack)  =  1.4  deg 
aej  (trim  elevator  deflection)  =  -2.36  deg 
6fj  (trim  flaperon  deflection)  =  -2.00  deg 
oTj  (nominal  engine  thrust)  =  5650,445  lb 


Longitudinal 

Body  Axis  Primed  Dimensional 

Derivatives 

Xq  =  -32.1643219 

4  = 

-.0008136 

M^  =  .0002939 

Xu  =  -.0157924 

Zu  = 

-.0000361 

Mu  =  -.0005463 

X,;  =  44.494278 

z:  = 

-2.0889177 

m;,  =  5.5969896 

Xq  =  -23.7776337 

.9999598 

Mq  =  -1.0057726 

X,',  =  -.6077153 
e 

-.2098653 

M',  =  -31.939163 
e 

X'.  =  19.4285533 

Of 

-.3693079 

M,',  =  -9.9644833 

Difference  Equation  Model  Parameters 

^1  = 

-3.9697145  0  i 

1 

Bi  = 

: .00206579 

1 

.00365134, 

0  -3.96971451 

I-. 3178785 

-.0992575' 

A2  = 

5.90880295  0  | 

B^  = 

i -.0062434 

-. 0109394 1 

0  5.9GSSC295; 

1.94639976 

.295502381 

A3  = 

i -3.90846236  0  | 

B3  = 

' .00622259 

.01090232; 

i 

!  0  -3.90346236' 

-.9401651 

-.2932328: 

A4  = 

1.9693739519  0  i 

' 

Bq  = 

!-. 002045 

-.0036142 

1 

1  0  .9693739519' 

1.31114387 

.09698786, 
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Table  A-3 


Longitudinal  Body  Axis  Primed  Dimensional  Derivatives 


=  ■ 

-31.8838959 

= 

-.0127706 

M,  = 

.0007531 

Xu  = 

-.0122359 

Zu  = 

-.0002911 

Mu  = 

.00006 

x«  =■- 

17.7885132 

z«  = 

-.4893466 

M.;  = 

1.7894831 

Xq  = 

-45.3607635 

Zq  = 

.9999136 

Mq  = 

-.3870971 

Xa  = 

1.6643476 

Z<3  = 

-.0770782 

m;  = 

-3.2519884 

e 

e 

e 

X')  = 

-4.29844 

^^f= 

-.0691358 

K  = 

Of 

.3253066 

Difference  Equation  Model  Parameters 

A 

1-3.99131 

1  ^  1 

0  1 
j 

Bi  = 

1.000768645  . 

j 

00068963' 

1 

1  0 

-3.99131  1 

i-. 03246486  . 

00324069! 

[5.97377273 

0  1 

62  = 

[-.00230457  - 

.0020662: 

A 

2  = 

i 

1 

1 

0  5 

.97377273: 

1 .097218426  - 

.0097182: 

[-3.9736152 

0 

63  = 

[ .002301487 

.0020636' 

3  = 

1 

0 

3.97361521 

1-. 09704231 

.0097143; 

[.991152575 

0  ! 

Bq  = 

I-. 00076556  - 

■.0006870' 

A 

4  = 

1 

1 

1 

0 

991152575: 

! .032288744  - 

•.0032363: 

Appendix  B 

Eguivalence  Bet'.veen  the  Step-Response  Matrix  and  the  Differepce 
Equation  N'odel's  Matrix  Coefficient  B| 


In  chapter  2  of  this  thesis  we  saw  how  the  use  of  the  step-response 
matrix  (  H(T)  )  provides  an  alternate  way  of  defining  the  control  law 
gain  matrices.  The  use  of  H(T)  a^so  provides  a  convenient  way  of 
relating  these  gain  matrices  to  the  matrix  coefficient  of  the  vector 
difference  equation  used  to  represent  the  aircraft  open-loop  dynamics. 
In  order  to  ilustrate  this  relationship,  it  is  convenient  to  consider 
the  significance  of  H(T)  and  how  it  relates  to  the  design  of  the 
fixed-gain  control  law. 

Consiu  first  the  solution  of  the  matrix  differential  equation 
given  in  Eq  (2-1),  and  the  output  relationship  given  by  Eq  (2-9).  For 
any  arbitrary  input,  the  output  response  is  given  by  (14:98-102) 


y(t)  =  C  exp(A  t)  x(0)  +  i  C  exp(A  (t-r))  B  u{-)  dr  (B-1) 

'0 


Under  the  assumption  of  zero  initial  conditions  for  the  plant  (x(0)  ~  0, 
y(t  t  0)  =0),  and  with  the  input  held  constant  between  samples  by  the 
controller,  the  response  of  the  plant  at  t  =  T,  for  a  given  u{0)  is 


y(T)  =  C  exp(A  (T-r))  B  t;(0)  d- 
^0 


(B-2a; 


=  H(T)  tv(O) 


(B-2b) 


Thus  H(T)  is  the  matrix  multiplier  that  relates  the  inp^t  u[0)  to  the 
output  of  the  plant  at  t  =  T. 

The  open-loop  plant  may  be  represented  by  a  vector  difference 
equation  form  given  in  Eq  (3-9).  For  the  same  inputs  and  initial 
conditions  as  before,  and  neglecting  the  zero  mean  equation  error  term, 
Eq  (3-9)  gives  the  following  output  relationship  at  t  =  T 

y(T)  =  Bi  y(0)  (e-3) 

If  the  output  in  Eq  (B-3)  is  to  be  the  same  as  that  of  Eq  (B-2),  for  the 
same  arbitrary  input,  we  must  then  conclude  that  H(T)  = 

The  question  of  how  H(T)  is  related  to  the  design  of  the  control 
law  gain  matrices  can  be  seen  simply  by  assuming  that  B  remains 
constant  between  any  to  given  sampling  instants,  and  doing  a  series 
expansion  on  the  matrix  exponential  in  Eq  (B-2a)  as  follows 

J 

y(T)  =  C  I  exp(A  (T-r))  dr  B  u(0)  (B-3a) 

^0 

=  C  f  {  I  +  +  H.O.T.  }  dr  B  u{0)  (B-3b) 

.Q  i.  z. 

For  small  sampling  periods,  the  contribution  of  the  terms  beyond  I  in 
the  series  is  relatively  small  and  may  be  neglected.  This  lead  to  the 
folowing  approximation  for  the  output  y(t) 

J 

y(T)  '  C  :  Id-  B  c'(0)  (B-4a) 

'0 

y(T)  5  T  C  B  u{0)  (B-4b) 
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This  apr-''Tix  is  ir,tsr',:ied  to  suple-'.ent  the  material  presented  in 
chapter  3  of  this  thesis  by  providing  additional  details  on  the 
implementation  of  tne  recursive  identification  algorithm.  The  first 
section  will  address  the  structure  of  the  difference  equation  model, 
while  section  two  will  emphasize  the  actual  implementation  of  various 
equations  in  the  algorithm. 

A.  St^'cture  of  t''e  Cn>~ference  Equation  Model 

As  described  in  Chapter  3,  the  difference  equation  model  used 
in  this  thesis  to  represent  the  plant  dynamics  has  its  origin  in  the 
transfer  ratr'ix  obtained  from  the  discretized  state  space  model  of  the 
plant.  The  transfer  natri/.  representation  leads  to  the  model  structure 
given  in  Eqn  (3-9).  In  ci'der  to  use  the  difference  equation  in  the 
estimation  algorithm,  it  miust  be  expressed  in  the  form  of  Eqn  (3-10). 
This  is  done  by  combining  all  the  delayed  input-output  measurements  in 
the  matrix  T'''(kT),  and  the  elements  of  the  parameter  matrices  A-j  and 
Ei  (i  =  1,2,..,N)  in  the  param.eter  vector  e.  In  forming  the 
parameter  vector  o  we  take  advantage  of  the  fact  that  the  matrices  A-j 
are  diagonal,  and  are  of  the  form  a^*!  (see  Appendix  A),  to  reduce  the 
number  of  elements  needed  in  c.  The  reason  for  the  form  af  the 
matrices  Aj  is  that  the  numbers  aj  represent  the  coefficients  of  the 
characteristic  equaticn  of  the  discrete  transfer  matrix.  Eqn  (3-10) 
then  assumes  the  form  given  in  Eqn  (C-l),  which  ilustrates  the  case  fcr 


a  second  order  model,  and  where  the  sampling  time  reference  has  t 
dropped  for  notational  convenience. 


yi(k)  0  0  '  -yi(k-l) 

yzik)  '  :  0  0  ui(k-l)  iy2(k-l)  -y2(k-l) 


.  /  \ . 


’  '  1 '  - ' 
:  tic:, 
bic 


♦  - 
*  -  /' 


(k-1) 


(k-2)  ;aiO.I) 

ib2(i.:; 

b2(2.1) 

b2(2,2) 

lanTlVl) 


i 

i  ^2(k)  i 


(C-1) 


The  "blank"  section  of  ■|■'^(kT),  underscored  by  (k-Z),  simply  ir::n:a:es 
that  the  arrangement  of  input-output  data  underscored  by  (k-1;  is 
repeated  as  many  times  as  required  by  the  order  of  the  difference 
equation. 

To  perform  the  scalar  measurement  update  of  the  parameter  vecccr, 
the  contribution  of  each  output  is  added  one  at  a  time.  This  is  dene 
assuming  that  the  equation  ei'ror  terms  in  Eqn  (C-1)  are  statist'cally 
independent  from  each  other.  This  is  considered  a  reasonable  assu-'pticn 
since  the  control  law  used  produces  decoupled  output  responses.  In  that 
case,  each  output  prediction  and  residual  is  generated  by  taking  the 
corresponding  rew  of  i''^(kT)  multiplied  by  the  cur'rent  estimate  or'  the 
parameter  vector'  c.  The  scalar  upoate  of  o  is  then  performed  erery 
time  step  kT  by  iter'ating  threugn  the  estimation  algorithm  as  man_,' 
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estimates  and  parameter  covariar- 
transfcrmaticn  that  avoids  tne  si 
expressions  replace  Eqns  (3-51) 
aljm'ithm  and  are  given  as  follcv.s 


:e  (i9:££)  by  a  simple  algebraic 
ngjlarity  ccnditicn.  Tne  resulting 
and  (3-65)  in  the  identification 


S'i(kT)  =  8'i.i(kT)  +  Pi-l^kT)  ..(kT)  r^fkT) 


(C-2) 


PT(kT)  =  Pi-i(kT) 


Fi_i(kT)  i(kT)  Pi-:i<T) 


vii<T)  ^  ;i(kT)  I  (C-3) 


V 


The  uphdte  eqeaticns  just  describes  are  ns.v  in  a  fcr::!  ir.sre  as-, enable 
for  computer  in^plerr.entation  but  further  r,;dlf icaticn  is  still  required. 
Kalrr,an-fi  Iter-type  update  equations  like  Eqns  (C-2)  and  (C-3)  are  renoun 
for  having  bad  nuT.erical  properties,  mainly  because  the  covariance 
update  equation  may  contain  differences  between  two  almost  equal  terms. 
Computer  roundoff  may  then  deteriorate  the  estimation. 

It  is  often  more  convenient  to  employ  factorization  techniques  for 
the  update  of  the  covariance  matrix.  A  commonly  used  factorization  is 
known  as  "U-D  factorization"  in  which  the  covariance  matrix  is  expressed 

'  ■  ' 

%  as  follows 


P(kT)  =  U(kT)  D(kT)  u'''(kT) 


(C-6) 


where  U(kT)  is  an  upper  triangular  matrix  with  all  diagonal  elements 
equal  to  1,  and  D(kT)  is  a  diagonal  matrix.  The  updating  of  the 
covariance  matrix  is  done  by  updating  the  U  and  D  fa:^:''5  rather 
that  updating  the  F  matrix  direccly.  The  U-D  factorization  has  good 
numerical  stability  and  enhances  the  algorithm's  tolerance  of  roundoff 
errors  (10,28,30  ). 

Detailed  discussions  on  the  U-D  algcriinm  can  be  found  in  >'6fer'ences 
(10),  (28),  (30),  (41)  and  (42).  Here  we  shall  limit  c.rsel.es  to  list  a 
segment  of  the  FCRTFAh  code  used  to  impl--enp  tn-  para--rte’'  esti-ate. 


and  covai'iance  matrix  updates  with  tn 
weighted  Gram-Schmict  ortr.cganal  izatic 
the  contrituticn  of  the  fa^it  detectar 
The  listing  is  limited  to  tncse  S‘ 


C3  25  KhCUT=l,h:.’ 

!  PEPFGP.M  THE  UPDATE  FOR  EACH 
!  OUTPUT  MEASUREMENT 

CO  30  1=2, lo:.: 

!  I  END  =  TOTAL  #  OF  PARAMETERS 

DO  30 

UTRAh(I,0)  = 

U^M1(0,I)  !  OBTAIN  U^i-iCkT) 

DO  35  I  =  1,IE“.0 

F(I)  =  0.0 

:  E. A. DATE  £d:i  i5-6Ca)  in  (19:39) 

DO  35  0=1,1 

!  f  =  U■•j_]^(^T)  ■•!(kT) 

F(I)  =  F(I)  + 

UTRAN(I,J)*GRA:mD(J, KNOUT) 

FTPFl  =  0.0 

DO  40  I  =  l,IEh'D 

G(I)  =  D(I)*F 

(I)  !  Eqn  (5-GCd)  in  (19:39) 

FTPFl  =  FTPFl 

+  F(I)-G(I)  :  T.(kT)  Pi-i(kT)  : i(kT) 

m 


C[-T~:,3-,LT;.<:,D.I^-.:.A,::,IlRR)  :  Eqn  (3-58) 

C-.L  A.~-:::,^Tr^:,X'.(vM.T),A8FnA,r,>.:jT),IERR)  I  Eqn  (3-62) 
ZETA  =  l.G  -  A.-r^;KMoT)*A.(r,M.T)  !  Eqn  (C-4) 


ETA,:)  =  a.\-:.:.T)  I  Eqn  (C-5).  Eqn  (5-6Cg)  in  (19:90) 
c:  A5  J  =  1,IE‘,0 

EETA(J)  =  BETA(J-l)  +  F( J)*G( J)*ZETA 


IF  (BETA(J)  .EQ.  0.0)  THEN  !  (19:365-366) 

!  Update  the  D(kT)  matrix 

D(J)  =  MIN(CREG,D(J)) 


D(J)  =  MIN(CREG,(BETA(J-1)/BETA(J))*D(J))  !  Eqn  (360h) 

!  in  (19:90) 

E'.OIF 


V.O)  =  G(J) 


IF  (BETA(J-l)  .EQ.  0.0)  THEN 
MU(J)  =  0.0 


!  Eqn  (5-E 


fTJ  ( J )  =  -F  ( J  ) '  T, 


4/4 


HD-MM  M7 
UHCUISSIFIED 


PMUMCTn-flOAPTIVE  HOOCL-FOLLONIM  FOt  IH-FLiaHT 
SIHUUITIOH<U)  flIR  FOtCC  INST  OF  TECH  MtlOHT-FRTTERSON 
NFt  OH  SCHOOL  OF  EHOIHCERIHO  L  R  FIHCIM  DEC  07 
RFIT/QE/EHB/07D-74  F/0  1/3.  3 


CONTINUE  !  with  next  measurement  update  ... 


TX’TV'.X.',  yC" 


*********************************************************************** 
*********************************************************************** 
SUBROUTINE  MWGS(B,N,UP,D,CREG,IERR) 

*****************************★****************★*■*■*■****■**■*■*:*■*****■*:****** 
*'*^:t'**7t'*********:i'*'**7t'5tTt*'*'7lr*'*'**'****5lr3t3lr*****5t5t'5k'******3t*X  Tt*5t'*******3lr****Tt'Tt*' 

DOUBLE  PRECISION  W(80,40) .UP(40.40) ,0(40) ,DB(80,TEMP , 

*  CREG,B(40) 

!  See  (28:333-334) 

DO  1  J=1,N 

DB(J)  =  D(J) 

DB(J+N)  =  1.0 
DO  1  1=1, N 

W(I.J)  =  UP(J.I) 

IF  (I  .EQ.  J)  THEN 

W(I+N,J)  =  SQRT(B(I)) 

ELSE 

W(I+N,J)  =  0.0 
qj',  END  IF 

1  CONTINUE 


DO  2  J=N 

,2,-1 

D(J) 

=  0.0 

DO  3 

K=1,2*N 

0(0) 

=  D(0)  +  W(K,0)*DB(K)*W(K,0) 

0(0) 

=  MIN(CREG,D(0)) 

DO  4 

1=1, 0-1 

TEMP 

=  0.0 

DO  5  K=1,2*N 

5  TEMP  =  TEMP  +  W(K,I)-^DB(K)*W(K,J) 


IF  (D(J)  .EQ.  0.0  THEN 


UP(I,J)  =  0.0 
GO  TO  4 


UP(I.J)  =  TEMP/D(J) 


ENDIF 


DO  6  K=1,2*N 


W(K.I)  =  W(K,I)  -  UP(I,J)*W(K,J) 
CONTINUE 


CONTINUE 


TEMP  =  0.0 


DO  7  K=1,2*N 


TEMP  =  TEMP  +  W(K,1)*DB(K)*W(K,1) 
D(l)  =  MIN(CREG,TEMP) 


RETURN 


-v*" v'"  v'  s’'*.'*- 


Root-Locus  Analysis  and  Time-response  Checks  with  Regards  to  Functional 

Controlabi 1 ity 

A.  Root-Locus  Analysis 

Early  in  the  design  stages  of  this  effort,  concerns  were 
brought  up  as  for  the  significance  and  effect  on  the  closed-loop  system 
behavior  of  a  transmission  zero  at  the  origin  that  leads  to  functional 
uncontrolability.  To  learn  more  about  this  condition,  a  root-locus 
analysis  of  the  of  the  closed-loop  system  was  conducted  with  the 
aircraft  model  corresponding  to  mach  0.9,  10,000  ft.  This  was  done  by 
individally  varying  the  design  parameters  a^,  c2,  and  p  in  the 
control  law,  and  observing  the  location  of  the  closed-loop  system  roots. 
To  simplify  the  problem,  actuator  and  sensor  dynamics  were  not  included 
in  the  analysis.  A  sample  of  the  results  is  shown  in  tables  D-1  through 
D-3.  In  these  tables,  R1  and  R2  are  the  fast  system  roots  and  R3 
through  R6  represent  the  slow  system  roots.  R5  and  R6  correspond  to  the 


Table  D-1 


Root-locus  Analysis,  gain  sweep  (c2  =  0.7  ,  p  =  0.8) 


roots 

ai  =  0.3 

o 

II 

i  n  =  0-5 

:  i 

Cl  =  0.7 

R1 

1 .297718585951 

'  .297714217964 

1 

1 .297705241223 

.286779451772 

R2 

; .637084434379 

.587157197317 

:  .437207541951 

.298180935417 

R3 

.991863119634 

'  .991392424014 

1  .991908999372 

.991930277332 

R4 

.9932CS391C24 

.993115733320 

' .993057794613 

.992988911001 

R5 

.999334953533 

.999534955393 

.999834957328 

.999334956957 

yw  MT  T*F^”f'T-^'\^'^r;‘.'; v/wv^Y,^»/  ». 


r-. 


*jr*  • 

■Cv 


X- 


Table  D-2 

Root-locus  Analysis,  02  gain  sweep  (ci  =  0.3  ,  p  =  0.8) 


roots^\ 

0 

II 

C\J 

(72  =  0*5 

(72  =  0.7 

02  =  0.8 

.596920647772 

.497341971943 

.297718585951 

.197827609203 

.687308842456 

.687155621738 

.687084484379 

.687070672164 

.991843993749 

.991853164897 

.991868119634 

.991874333508 

.993793593327 

.993522421431 

.993208391024 

.993108852556 

.999847457184 

.999841354479 

.999834953500 

.999833067057 

1.0 

1.0 

1.0 

1.0 

Table  0-3 

Root-locus  Analysis,  p  gain  sweep  (c^  =  0.3  ,  02  =  0.7) 


roots''^^ 

P  =  0.6 

p  =  0.7 

P  =  0.8 

II 

0 

R1 

.295710430846 

.296713062375 

.297718585951 

i  .298727026744 

R2 

.685130298700 

.686104281008 

.687084484379 

'  .688071029587 

R3 

.993924862324 

.992898588723 

.991868119634 

; .990833375662 

R4 

.995108546177 

.994161375687 

.993208391024 

j .992249858364 

R5 

.999840396440 

.999837226695 

.999834953500 

1.999833244131 

'  R6 

1 

1.0 

1.0 

1.0 

1  1.0 

roots  introduced  by  the  vector  integrator  of  the  control  law. 

As  shown  in  these  tables,  the  roots  R1  and  R2  approach  the 
asymptotic  values  defined  by  Eqns  (2-41)  and  (4-9)  given  by 


\  i  =  1  -  c  ■ 
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The  slow  roots  R3  and  R4  approach  those  of  Eqn  (2-39)  which  by  virtue  of 
Eqn  (2-44)  reduces  to 

Xt  =  1  -  T  p  i  =  3,4  (D-2) 

while  the  remaining  slow  roots  R5  and  R6  correspond  to  the  roots  of  Eqn 
(2-40).  In  this  case,  R5  approaches  a  transmission  zero  of  the  plant 
located  at  0.99982119878,  while  R6  sits  on  top  of  the  transmission  zero 
at  the  origin  produced  by  the  selection  of  pitch  rate  as  an  output 
quantity.  It  is  obvious  that  the  pole-zero  cancelation  occuring  between 
the  transmission  zero  at  the  origin  and  a  pole  from  the  vector 

integrator  does  not  allow  this  system  root  to  be  placed  arbitrarily. 

This  situation  renders  the  system  functionally  uncontrollable,  however, 
this  is  to  be  expected  given  the  kinematic  definition  used  in  the 
state-space  model  of  the  aircraft.  The  status  of  functional 

uncontrolabi  1  ity  in  this  example  does  not  necessarily  implies  that 
satisfactory  tracking  performance  may  not  be  achieved,  but  simply 
expresses  the  invariability  of  the  kinematic  relationship  between  the 
pitch  angle  and  pitch  rate  responses.  Throughout  the  analysis,  all 
system  roots  behaved  in  a  predictable  way  and  remained  at  or  within  the 
unit  circle  as  desired.  Satisfactory  tracking  performance  can  then  be 
achieved  provided  one  considers  the  kinematic  relationship  just 
mentioned  when  generating  the  input  commands  to  the  closed-loop  system. 

B .  Time  Response  Checks 

Another  concern  regarding  the  issue  of  functional 

controlabi 1 ity  was  wether  cr  not  the  closed-loop  system  was  capable  of 


Since  the  command  Inputs  used  throughout  this  thesis  were  not  the  result 
of  decoupled  motion  from  the  AFTI/F-16,  it  may  be  argued  that  these 
signal  are  "related"  to  each  other  and  thus  not  address  the  situation  in 
question.  To  that  effect,  the  command  inputs  were  modified  and  applied 
to  the  same  aircraft  model  of  part  A  with  the  nominal  control  law  design 
of  chapter  4.  The  results  are  given  in  Figures  D-1  through  D-32. 

Figures  D-1  through  D-8  present  the  case  where  the  flight  path  angle 
command  is  temporarily  "clipped"  to  hold  a  maximun  value  of  2  deg.  while 
the  pitch  rate  command  remains  unchanged.  The  responses  show  no 
noticeable  tracking  performance  degradation. 

Figures  D-9  through  D-16  show  the  responses  to  a  fligtit  path  command 
of  0  deg.  while  the  pitch  rate  command  again  remains  unchanged.  In 
this  example,  the  control  law  contained  no  anti-windup  provisions. 
Figures  D-9  and  D-10  show  relatively  good  tracking  performance  for  the 
flight  path  angle  response.  However,  the  pitch  rate  response  suffers 
from  a  short  period  of  instability  caused  by  deflecting  the  flaperon  to 
its  position  limit.  The  larger  control  surface  deflections  are 
attributed  to  the  bigger  demands  on  control  action  imposed  by  the 
decoupled  maneuver.  Aside  from  the  instances  where  the  flaperon  is  at 
its  deflection  limit,  the  responses  exhibit  good  tracking  behavior. 
Thus  the  condition  of  functional  uncontrolab i  1  i ty  is  not  seen  as  the 
cause  of  the  problem. 

The  previous  maneuver  was  repeated,  this  time  with  the  integrator 


limiter  mentioned  in  chapter  4.  This  is  to  test  the  effectiveness  of 
the  antiwindup  protection  scheme.  The  results  are  shc^n  in  Figures  D-17 
through  D-24.  These  ilustrate  that,  althougn  the  flaperon  deflection 


Figure  D-1.  Flight  path  angle  command  and  response  (deg). 
"Clipped"  flight  path  command. 


Figure  D-9.  Flight  path  angle  command  and  response  (deg). 
Flight  path  command  =  0  deg.  No  integrator  limiter  in  control  law 


Figure  D-10.  Flignt  patn  angle  ti'acking  performance  index  (deg). 
Flight  path  cc—and  =  0  deg.  No  integrator  limiter  in  control  law 


Figure  D-13.  Elevator  deflection  (deg). 

Flight  path  command  =  0  deg.  No  integrator  limiter  in  control  law 


Figure  D-14,  Ele'.ator  deflection  rate  (deg/sec). 

Flight  path  command  =  0  dec.  No  integrator  limiter  in  control  law 


T  :r-e  (  ££c  ; 


Figure  D-19,  Pitch  rate  ccT-and  and  response  (deg/sec). 
Flight  path  corr-nad  =  0  deg.  Control  la.v  with  integrator  limiter 
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Figure  D-21.  Elevator  deflection 
Flight  path  cor’.xard  =  0  ceg.  Control  law  with 


(deg). 

integrator  limiter 


{  ceg  'sec ) . 
inCtr^ratcr  limiter 


reaches  its  limit,  in  this  cccasicn  the  amc^nt  of  time  that  the  surface 
remains  at  its  rraomjm  deflecticn  is  ccnsicerably  less  than  in  the 


pr'evicus  simjl; 


Ccr.setuently ,  tr.e  aircraft  repcnses  exhibit 


imprcved  behavior  with  regards  to  tne  shorter  duration  of  the 
instability,  and  tne  g:;d  tracking  perfc'-mance  that  follows  once  the 
cc"'"'and-d  control  surface  deflections  fall  tack  within  the  given 
constraints  for  this  example. 

Finally,  the  system  is  tested  with  the  pitch  rate  command  set  to  0 
deg/sec,  while  the  flight  path  command  is  set  as  in  the  original 
maneuver.  The  integrator  limiter  is  also  used  in  this  simulation.  The 
aircraft  responses  fer  this  simulation  are  are  presented  in  Figures  D-25 


hreuen  D- 


Tracking  performance  is  considered  good,  despite  the 


P'"esen:e  of  a  s^all  transient  in  the  pitch  rate  response  caused  by  a 
brie'  saturation  of  tne  flapercn  deflecticn. 

C.---’'a'l,  these  ti"'e  respc'ses  indicate  trat  gi.en  sufficient  control 
p:.’.-;’'  and  cc"-an;d  ra",':„vei's  trat  are  witnin  t‘ie  capacilities  of  the  host 
vehicle,  the  issue  cf  functicral  unccntrclaci 1 ity  {as  it  applies  in  this 
f'-rsis)  dees  net  precludes  tne  tracking  cf  ceccupled  com.mand  inputs. 
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In-flight  simulations  are  normally  accomplished  by  using 
model-following  control  laws  which  depend  on  accurate  knowledge  of  the 
stability  derivatives  of  the  host  aircraft.  Degraded  simulation 
performance  may  result  if  the  stability  derivatives  deviate  considerably 
from  their  presumed  values.  Gain  scheduling  is  often  employed  to 
compensate  for  plant  parameter  variations,  but  this  form  of  open-loop 
compensation  usually  requires  extensive  flight  testing  for  proper  fine 
tuning. 

This  thesis  implements  an  adaptive,  fast-sampling  control  law  to 
compensate  for  changing  aircraft  parameters.  The  step-response  matrix 
which  is  required  for  this  implementation  is  identified  recursively 
using  a  recently  developed  technique  which  does  not  require  special 
"test"  signals  and  which  automatically  discounts  old  data  depending  on 
the  input  excitation  detected.  Tracking  fidelity  is  maintained  despite 
parameter  changes  which  occur  either  abruptly  or  slowly.  Simulations 
are  conducted,  using  a  model  of  the  AFTI/F-16  aircraft  and  the  control 
design  package  MATRIX^,  to  test  the  resulting  adaptive  system.  Actuator 
position  and  rate  limits  are  discussed.  The  performance  of  the 
resulting  system  is  excellent  and  demonstrates  the  relative  advantages 
of  adaptive  controllers  for  in-flight  simulation.  Recommendations  are 
made  for  future  analysis  including  the  use  of  moving-bank  estimators. 


